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Abstract 

The behavior of a Lattice Monte Carlo algorithm (if it is designed correctly) must approach that 
of the continuum system that it is designed to simulate as the time step and the mesh step tend 
to zero. However, we show for an algorithm for unbiased particle diffusion that if one of these 
two parameters remains fixed, the accuracy of the algorithm is optimal for a certain finite value 
of the other parameter. In one dimension, the optimal algorithm with moves to the two nearest 
neighbor sites reproduces the correct second and fourth moments (and minimizes the error for the 
higher moments at large times) of the particle distribution and preserves the first two moments 
of the first-passage time distributions. In two and three dimensions, the same level of accuracy 
requires simultaneous moves along two axes ("diagonal" moves). Such moves attempting to cross 
an impenetrable boundary should be projected along the boundary, rather than simply rejected. 
We also treat the case of absorbing boundaries. 
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I. INTRODUCTION 



Particles diffusing in a liquid medium perform random walks. In a computer simulation 
of this process, it is sometimes convenient to introduce a lattice, assume that the particles 
can only reside at lattice sites, and allow only certain types of moves between the sites 



171 ]. Different variants of such Lattice Monte Carlo (LMC) algorithms are distinguished by 
the sets of allowed moves and the corresponding probabilities. These are chosen depending 
on the purpose of the algorithm. For instance, Metropolis Monte Carlo lSj can be used 



to study equilibrium properties, but is in general inadequate for dynamics, especially in 



a strong external field 12j. For dynamical algorithms, another important (and often ne- 



glected) ingredient is the time step (the amount by which the time is incremented after each 
attempted LMC move) [19j. While in principle the time step can vary during a simulation, 
in this paper we restrict ourselves to algorithms where it remains constant. Moreover, we 
only consider the case of unbiased diffusion, where there is no external field or other exter- 
nal influence (such as a flow) that would drive the particles preferentially in a particular 
direction. Some of the approaches used in this paper are applicable to the case of biased 
diffusion as well, and this will be described in a separate publication. We also assume that 
the system is uniform, i.e., the diffusion constant is the same everywhere. 

As an example of a dynamical LMC algorithm, consider unbiased one-dimensional (ID) 
diffusion along an infinite line, with the ID lattice sites placed equidistantly and numbered 
consecutively by integer numbers. In the simplest approach, at each step the particle moves 
left or right to a neighboring site with equal probabilities. If for a particle that is known to 
be at site i, the probability that it moves to site j during a given time step is p^j, then 

/', ■:, ■ ■ i>j -j i = V 2 , (1) 
= Ofor l^j±l. (2) 

The time step r can be determined, e.g., using a mapping onto the mean first-passage time 



between sites [19], and is equal 



a 2 , . 

where a is the distance between the lattice sites (the lattice constant or the mesh step) and 
D is the diffusion constant of the particle in the medium. 

Rather than working with individual particles, one can look at the evolution of particle 



distributions, described by the set of the mean particle numbers at each site, {rii(t)}, where 
the subscript i refers to the site number. In general, the mean particle numbers after the 
move are given in terms of those before the move by 



(*+T) = Z)p'-tf n '(*)' ( 4 ) 



where in principle 1 1 — j\ can be allowed to be larger than unity. Equation (prj is known as 
the master equation. For the particular simple algorithm described by Eqs. ([I])-©, 

n j {t + r) = -{n j - 1 {t)+n j+1 {t)). (5) 
The set of equations (jSJ) for all sites can be solved numerically for given initial conditions, 



and this provides, in a sense, a numerically exact solution of the LMC algorithm 20]. 

On the other hand, diffusion can be studied using continuum equations that describe the 
time evolution of the particle concentration. For example, for unbiased diffusion in ID 

dn(x, t) = d 2 n(x,t) 

dt ~ dx 2 ' [ ' 

To solve such an equation numerically, one can discretize it replacing the derivatives with 
differences: 

n(x,t + r) — n(x,t) n(x + a, t) — 2n(x, t) + n(x — a, t) , . 

~ D x - (7) 

r a 1 

or 

Dt 

n(x, t + t) « n(x, t) H — (n(x + a,t) — 2n(x, t) + n(x — a, t)) . (8) 

If x = ja, where j is integer, then we get a set of equations involving only the values of 
n(x, t) at discrete points: 

Dt 

nj (t +t) w nj (t) + — (n j+1 (t) - 2 nj {t) + %-!(*)) , (9) 

(X 

where nj(t) = n(ja,t). Note that this coincides with Eq. (jSJ) when r is given by Eq. ([3]). In 
other words, the master equation for a LMC algorithm can be viewed as a discretization of 
the continuum equation for the same diffusion problem. 

Note, however, that other choices of the time step r are possible in Eq. (Q. For each 
such choice, one can define a new LMC algorithm with a new set of transition probabilities 
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FIG. 1: A schematic of the ID algorithm considered in this paper. At a particular step, in addition 
to the usual moves to nearest neighbors, the particle can also stay put (which can be considered 
as the move to the original site). The notation for the probabilities of the moves is given. These 
probabilities (along with the time step) can be adjusted to optimize the algorithm. For unbiased 
diffusion, we always have pj_i_^- = Pj+i^j- 



by comparing Eq. to Eq. (j4]) (with the caveat that these probabilities have to remain 
bounded between zero and one). Indeed, if we rewrite Eq. (Q as 

x Dt , s / 2Dt\ . . Dt 
n 3 {t + r) « — nj-iit) +1 r )n j (t) + —n 3+l {t) (10) 

and compare this expression to Eq. PJ, we obtain 

Pj : -, /'/ •: •, = "^"> I 11 ) 

2Dt . . 

Pj-+i = l -^T- ( 12 ) 

To ensure that these probabilities remain between zero and one, the time step must be in 
the range 

0<r<^. (13) 

Note that in general we now have a nonzero probability Pj->j of staying at the same site 
during a particular time step (Fig. [p. The only exception is the particular case of the 
algorithm given by Eqs. (JTJ) — (J3J) , which corresponds to the largest possible time step, 

W = (M) 

Such an algorithm (to which we will refer as the ordinary algorithm) is a logical choice from 
the point of view of the computational efficiency of the simulation. However, choosing a 
shorter time step can increase the accuracy of the LMC algorithm. 
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The problem of optimizing the accuracy of the algorithm given by Eqs. (TIT]) and (fl2|) seems 
trivial at first. Indeed, the smaller the time step and the mesh step, the more accurately 
the discretization represents the continuum equation §6§ and thus the more accurate the 
corresponding LMC algorithm. In practice, however, computational resources are limited 
and one can only choose the smallest time step for which the simulation is still feasible. 
The practically relevant question then is: given that time step, what is the optimal mesh 
step? Based on Eq. ffl3|) . it cannot be infinitely small. The smallest possible mesh step (the 
finest space discretization), a min = a/2Dt, corresponds to the ordinary algorithm described 
by Eqs. (|TJ) (J3j) . While at first glance it seems it should provide the best accuracy, we will 
see that this is not the case. Conversely, in some cases one may wish to fix the mesh step 
(for example, choosing the largest value that ensures that spatial inhomogeneities are still 
reproduced accurately) and then optimize the time step. The largest possible time step, 
t = a 2 /2D [again, corresponding to Eqs. ([1]) and (|2J)], as mentioned, is optimal in terms of 
computational speed. As for the accuracy, at first glance it may seem that the smaller r the 
better, but again, this trivial answer turns out to be incorrect. 

This problem of optimizing the accuracy of LMC algorithms for diffusion given a fixed 
time step or mesh step is the subject of this paper. While in our considerations we fix the 
mesh step and allow the time step and the transition probabilities to vary, it will be clear 
from our derivation that fixing the time step instead will lead to the same probabilities and 
the same relation between the mesh step and the time step at optimality, at least as long 
as all inhomogeneities, such as obstacles and their features, are much larger than the mesh 
step. In the next section, we solve the optimization problem for the ID diffusion algorithm 
described above. We do this using three different approaches, with the same end result. 
In Sees. IIHI and IIVt we consider the 2D and 3D analogs of the same problem and show 
that the same accuracy as in ID requires the introduction of simultaneous moves in at least 
two directions, or diagonal moves. In Sees. |V] and I VI I we show how the LMC rules should 
be modified for sites next to reflecting or absorbing boundaries. We end the paper with a 
discussion of our results. 
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II. ONE DIMENSION 



In this section, we optimize the ID LMC algorithm for unbiased diffusion with the moves 
restricted to nearest neighbors. We do this in three different ways. The first approach is 
based on the comparison between the solutions of the continuum equation (E]) and the master 
equation (j3j) (with only £>j±i_>j and Pj->j nonzero); they should match particularly for long- 
wavelength, slowly decaying modes. The other two methods compare the exact moments 
of the particle distribution or the first-passage time (FPT) distribution with those of the 
respective distributions generated by the algorithm. Using these three different approaches 
reveals more properties of the optimal algorithm, but the resulting sets of parameters are the 
same. We also consider the full FPT distributions and show numerically that while for both 
the ordinary algorithm [Eqs. (JTJ) — (J3J)] and the optimally accurate algorithm the distribution 
eventually approaches the correct one as the mesh gets finer, the latter algorithm achieves 
the same accuracy for a much coarser mesh. 

A. Comparison to the continuum equation 

Note first that the general solution of the continuum equation ([6]) on an infinite line can 
be written as a sum (or, rather, an integral) over modes with different wave numbers k, each 
of which decays exponentially in time: 



where C(k) is an arbitrary complex function [except for the fact that we need C(—k) = C*(k) 
for the solution to be real] and the dispersion relation is 



The general form of the master equation for a LMC algorithm in ID where only moves 
between neighboring sites are allowed is 



The general solution of the system of master equations (fTTj) can be written in a form similar to 
Eq. (|T5|) . except for the fact that the integration limits are finite due to space discretization: 




(15) 



a c {k) = Dk 2 . 



(16) 



7ij(< + r) = pj^jUjit) +pj_i_^n,-_i(t) +Pj+i^jn j+1 (t). 



(17) 




(18) 
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The dispersion relation ad(k) can be obtained by substituting Eq. (1T81) in Eq. (1TT1) : 

ad{k) — In (pj^j + Pj-i->j exp(-ika) + Pj+i->j exp(ika)) . (19) 

T 

Note that Eq. (fT9~|) is the esaci dispersion relation for the discrete equation ffT7|) . even for 
finite values of a and r. 

The closer Eq. (119]) approximates the dispersion relation (116]) for the exact, continuum 
equation, the more accurate the discretization and the corresponding LMC scheme. One 
can argue that the behavior of the solution of the diffusion equation on the longest length 
scales (smallest k) should be reproduced in the first place, especially given that these longest- 
wavelength modes also take the most time to decay. For this reason, we expand the dispersion 
relations [Eqs. (116]) and (119]) ] in powers of k and try to match as many terms as possible. 
Since Eq. (116]) contains only a k 2 term, we need to find the conditions that would eliminate 
as many terms as possible in the series expansion of Eq. (119]) . starting with the lowest order, 
except for the k 2 term whose coefficient should equal D. The k° term disappears if 

Pj^j + Pj-i^j + Pj+i^j = L (2°) 

In other words, the probabilities must be normalized. Given Eq. (1201) . we find that the k 1 
term in Eq. (119]) has a zero coefficient if 

— (Pi-i->j -Pj+i^j) = 0- (21) 

This means that the probabilities of moving to the left and right should be equal, which 
makes sense since our diffusion process is unbiased. Taking conditions (1201) and (12T1) into 
account, Eq. (1T91) becomes 

a d {k) = — In (1 + Pj _ !^A-k 2 a 2 + fc 4 a 4 /12 - fc 6 a 6 /360 + ...))■ (22) 

T 

After expanding the logarithm in a series, the k 2 term of the resulting expression matches 
Eq. (USD if 

= D . (23) 

r 

This is equivalent to Eq. ( ITT]) and ensures the correct diffusion rate. Requiring that the next 
(A; 4 ) term be equal to zero gives the expression 

_ip^ + M!) =0 . (24) 



The solution of this equation is simply 

Pi-i-tf = g- ( 25 ) 

Note that pj-i^j = is not a solution, since in the limit Pj-i^j — > not just the expression 
in the parentheses in Eq. f l24|) . but also r approaches zero and the ratio Vi-\^il T remains 
finite, according to Eq. (|23l) . This means that the k 4 term is not eliminated in the limit r — > 
and so, perhaps surprisingly, the algorithm is not optimal in this limit. From Eqs. ( |20l) . ( I2TT) . 
and (ESD, we then get 

Pi±i-+i = \, (26) 

Pi-tf = | ( 2 7) 
a 2 . , 

T = 55' (28) 

Interestingly, the optimal algorithm (in terms of accuracy) requires a time step that is 1/3 
of the maximum value r max allowed for a particular step length a (as given by Eq. (!T4l)); as 
a consequence, the particle must stay put 2/3 of the time. While this is costly in terms of 
computing time, it is obviously better than the naive expectation that optimality would be 
achieved for r — > even when a is finite. 

Ideally, further terms in the series expansion of Eq. ffl9|) should also be zero. All odd- 
order terms are automatically zero, but for even-order terms starting with 0(k G ) this is 
impossible to achieve since we have run out of adjustable LMC parameters. Indeed, using 
Eq. ( 1221) with Pj-i^j = 1/6, one can check that the k & term is nonzero. Of course, this 
term (and all subsequent terms) do vanish in the limit a — > (which also automatically 
means r — > 0). For a fixed and finite mesh step a, one can introduce more parameters by 
allowing longer-range moves, but this can cause problems since it would effectively increase 
the coarseness of the discretization of space (e.g., with jumps of size 2a, a particle could 
move through an obstacle of size a). 

In Fig. [2J we compare the exact dispersion relation a c (k), as given by Eq. ( fl6|) . with the 
dispersion relations that correspond to LMC algorithms built with different time steps r. 
The latter dispersion relations are obtained using Eq. (fl9|) with the probabilities given by 
Eqs. fTTTl) and (Tl2l and thus satisfying Eqs. fTSOj) . fT2Tl) and (T231 . which guarantees the correct 
expansion terms up to k 2 , but the fc 4 term only vanishes for r = |r max . For the maximum 



S 



possible value r = r max = a 2 /2D, which corresponds to the ordinary algorithm with zero 
probability of staying on the site, two peculiarities catch the eye. First, at k — 7i/2a, ctd(k) 
diverges, i.e., the corresponding mode decays infinitely fast. Indeed, the distribution 

...,0,1, 1,1, 0,1 1,1 0,... (29) 

decays to the uniform distribution . . . , |, ~, |, . . . in a single time step and does not change 
afterwards. Even more strikingly, for k = ir/a, the real part of a<i(k) is zero, i.e., this mode 
does not decay at all; at the same time, the imaginary part is Im ad^/a) = n/r. As a 
result, the initial distribution 

...,0,1,0,1,... (30) 

oscillates indefinitely: 

...,0,1,0,1,... 1,0, 1,0,... -K..,0,1,0,1,...-K.. (31) 

One practical consequence of this is that a particle starting at one of the even-numbered 
sites will always visit only odd-numbered sites at odd-numbered time steps, and only even- 
numbered sites at even-numbered time steps. When a smaller time step is used with the 
same lattice constant a, there is some probability of staying on the same site at each time 
step, and oscillations similar to Eq. ( I3T1) decay. Eventually, for r < |r max , the oscillations 
do not occur at all, as ad is always real in that case. The exact dispersion curve is matched 
best for small k when r = |r max , as expected. 

Instead of using the explicit solution of the master equations [Eq. (fl8|) ]. another possible 
approach is to insert the partial solutions of the continuum problem [i.e., the integrand 
of Eq. ( fT5l) ] in the master equation (TT7]) replacing x with ja and find for what values of 
the parameters the resulting equality is satisfied most accurately for small k. After the 
substitution, dividing both sides by C(k) exp(ikaj — Dk 2 t), we obtain 

? 

exp(—Dk r) ~ pj^j + Pj-\-+j exp(-ika) + Pj+i^j exp(ika). (32) 

Here we use the ~ sign as a reminder that this is a relation that needs to be verified (hence 
the question mark); moreover, it is not expected to hold exactly, since lattice diffusion is only 
an approximate representation of the continuum diffusion process (thus "~"). By expanding 
both sides in the Taylor series in k we can verify that this equality is satisfied to 0(k 4 ) only 
when the probabilities and the time step are given by Eqs. (T26]) — (T28]) . 
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FIG. 2: For Fourier components of a one-dimensional particle distribution, characterized by their 
wavenumber k, a comparison between the exact decay rate obtained by solving the continuum 
diffusion equation [given by Eq. f)16|> : thick line] and the decay rate observed in a LMC algorithm 
with time step r [given by Eq. (|19j> ] . for several different values of r (thin lines). For the maximum 
possible time step, r = r max = a 2 /2D, the decay rate diverges at ka = ir/2, and for larger ka it 
is complex becoming purely imaginary (Re a = 0) at ka = it. As t decreases, the wavenumber 
at which the divergence occurs shifts to higher k reaching ka = it when r = r max /2. The closest 
matching between the exact and LMC decay rates at small k is observed when r = r max /3. 



An algorithm with a nonzero probability to stay put was proposed before 19] for biased 
diffusion in an external field. It was shown that in that case this nonzero probability was 
necessary even to reproduce the correct diffusion constant. Interestingly, the zero- field limit 
of that algorithm coincides with the optimal algorithm derived here, even though in zero 
field even the ordinary algorithm reproduces the diffusion constant correctly. 

B. Moments of the particle distribution 

Another way of constructing and analyzing LMC algorithms is by looking at the moments 
of the particle distribution. In continuum space, for particles starting at the origin (x = 0) 
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at t = 0, the distribution at time t is 



n(x,t) = -^= e - x2 / 4Dt . (33) 

The moments of this distribution are 

(x 2m ) = (2m - l)\\(2Dt) m = (2m - 1)!! x (x 2 ) m . (34) 

All odd moments are zero. The first three nonzero moments are 

(x 2 ) = 2Dt, (35) 
(x 4 ) = 12DH 2 , (36) 
(x 6 ) = 120D 3 t 3 . (37) 

We will now compute the moments of the particle distribution during a LMC random 
walk. For a lattice walk starting at site 0, the position after N steps is 

N 

x N = a^r/i, (38) 

i=l 

where rji is +1 (—1) for a move to the right (left) and when the particle does not move. 
We assume from the outset that Pj-i-^j = Pj + \^j, which makes all odd moments zero 
automatically. Before proceeding, it is convenient to note that for any positive integer m, 

( V 2m ) = 2m x Pj ^ + l 2m x Pj _^ + (-l) 2m x p j+1 ^ 

= I yj • Pj • I >j = 2pj-l->j, (39) 

(^-i) = O 2 " 1 - 1 x Pj ^ + I 2 ™" 1 x p 3 _^ 3 + (-I) 2 ™" 1 x p j+1 ^ 

= Pj-l-tf ~ Pj+l->3 = 0- ( 4 °) 



The second moment, which is the average square of Eq. (138)) . is 



N 



\ X N 



) = a 2 J2^Vk). (41) 



i,k=l 

Since different steps are uncorrelated, 



(ViVk) = S ik (i] 2 ) = 25 ik Pj-i-*j. (42) 
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Here and below 8^...^ is (the generalization of) the Kronecker's delta, which is unity when 
all indices coincide and zero otherwise. We obtain 

X - ,2. 



( X 



2a v 

2a 2 p J ^ J § ik = 2a 2 p J . 1 ^ j N = ^±^ t. (43) 



r 

i,fc=l 



This coincides with the continuum result, Eq. ( 1351) . when a 2 pj-i^j/r = D. This is the same 
as Eq. (TTTT) [or (|23|) ] that is obtainable from matching the k 2 terms in the continuum and 
discrete dispersion relations. Likewise, 

N 

(x%) = a i (ViVkVrVs)- (44) 

i,k,r,s=l 

The average (r]ir] k ri r ri s ) equals: 1) (rf) 2 when there are two pairs of equal indices, but the 
indices in different pairs are different; 2) (f? 4 ) when all four indices are equal; 3) otherwise. 
As a single expression, 

(ViVkVrVs) = (S ik S rs + 5 ir 5 ks + 5 is 5 kr - 35 ikrs )(r] 2 ) 2 + 5 ikrs (r] 4: ) 

= A(5 ik 5 rs + 5 ir 5 ks + 5 is 5 kr - 35 ikrs )p 2 _ 1 _ !fj + 28 ikrs pj^ 1 ^j 1 (45) 



and we get 



(x%) = a 4 [(l2N 2 -12N)p 2 j _ 1 ^ j + 2Np j ^ j ] 

An o/j-2 /_2\„2 , /i/_\/o„ ioJ 



= oWWi^ + (tM(2Pj-i-+ S - 12p^-)]- (46) 
When Eq. (TTTj) is satisfied, so the second moment is correct, this becomes 

( 4) = 12flV + C*-^ - ifr?-.~>y (47) 



This approaches the continuum result, Eq. fl36|) . when t — > oo, but only coincides with it 
exactly at all times, if Pj-i^j = 1/6, which coincides with Eq. (!25l) and gives rise to the 
optimal algorithm [Eqs. fl26|) - fl28|) ]. 

In general, even moments are given by 

N 

(x 2 N m )=a 2m J2( Vn ... Vl2m ). (48) 



i,=l 



For the sixth moment, we need (rj^ . . . f]i 6 ). This equals: 1) [f] 2 ) 3 when there are three pairs 
of equal indices, but no equality of indices between any of the pairs; 2) (rj 2 )(r] 4 ) when there 
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are four equal indices and another pair of equal indices not equal to the first four; 3) (rj 6 ) 
when all six indices are equal; 4) otherwise. When the indices run from 1 to N, there are 
15N(N — l)(iV — 2) combinations of indices of the first type, 15N(N — 1) combinations of 
the second type, and N combinations of the third type. Then 

( x e N ) = a e [15N(N- l)(iV-2)(77 2 ) 3 + 15N(N - l)^ 2 )^ 4 ) + iV(r/ 6 )] 

= a 6 [120N(N - 1)(JV - 2)p)_ 1 ^ + Q0N(N - l)p)^ + 2N Pj ^\ 
= a 6 [120(t/r) 3 ^ 3 _ 1 ^. + (t/r) 2 (60p 2 _ 1 ^. - 

+(t/r)(2p i _ 1 _> i - Q0p 2 _^ + 240pj_ 1 _+ j )). (49) 

It is easy to check that the coefficient of the t 3 term matches Eq. (157)) whenever Eq. ( II ip 
is satisfied, i.e., whenever the second moment (and the leading term in the fourth moment) 
are correct. Moreover, the t 2 term is zero if and only if Pj-i^j = 1/6, i.e., for the optimal 
algorithm. However, even in this case the remaining t 1 term is incorrect. Yet, the error 
of the sixth moment at large t is optimized, since in this case the relative error is 0(t~ 2 ), 
whereas in all other cases the t 2 term is present and the error is 0(t _1 ). 

In fact, this statement about the optimization of the error is true for all moments. In 
the case of the (2m)-th moment, the two terms of the highest order in t are of order t m and 
t m_1 . By analogy with the sixth moment, the only combinations of indices in (77^ . . -T]i 2m ) 
contributing to these terms are: 1) those with m pairs of equal indices with no equality 
between pairs; 2) those with four equal indices and m — 2 pairs of equal indices with no 
equality between these groups. For the first type, there are (2m — 1)!! ways to break the 
indices into pairs and N\f(N — m)\ ways to assign the values of the indices to these pairs, 
for the total of (2m — 1)\\N\/(N — m)\ terms in the sum, each of which is (rj 2 )'" 1 . For the 
second type, there are ( 2 ™) ways to choose the group of four indices, (2m — 5)!! ways to 
break the remaining indices into pairs, and N\/(N — m + 1)! ways to assign the values of 
the indices to all these groups, for the total of ( 2 ™)(2m — 5)!!iV!/(iV — m + 1)! terms, each 
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of which is (i] 4 ) (rf) m 2 . The (2m)-th moment then is 



( x 2 ™) = a 2m 



(2m-l)!!JV! 



,2m 



(N -m)\ 
(2m- 1)!! ( N 



{vT + 



! ™)(2m-5)!!iV! 
(N-m + 1)! 



(Vliv) +0(N m ~ 2 ) 



m(m - 1 V-M(2p,- 1 -,„ 



+ ^ 2 ^ (2m - 5)\\N m ' 1 (2p j _ 1 ^ j ) m ' 1 + 0(iV m - 2 ) 



,2m 



(2m- 1)!! 



2p, 



+ 



771 — 1 jn 777, 

(2m - l)!!m(m - ^^ hzlzliLZlizk^i ^-i + 0(r -2) 



(50) 



We see again that the coefficient of the leading t m term is correct whenever Eq. (fTTl) is 
satisfied, and the subleading term vanishes when Pj-i-^j = 1/6. 



Moments of the FPT distribution 



Consider again a particle diffusing in continuum space and starting at x = at time 
t = 0. In the FPT problem 2l| . we consider two imaginary walls at x = ±6 and find the 
first instance at which the particle reaches one of these walls. The corresponding time T is 
the FPT. The LMC analog of the FPT problem would be starting at site at t = and 
determining the first time one of the sites ±iV is reached. In a good LMC algorithm, the 
corresponding FPT distribution should match the continuum one for b = Na as closely as 
possible. 

Consider first the case N — 1 (and thus b = a). In the LMC, the FPT then corresponds 
to the step at which the particle first moves out of site (and into one of its neighboring 
sites ±1). This happens at mth step with probability 



Then the first moment of the FPT, or the mean FPT (MFPT), is 



(51) 



T 



m7i„ 



a 

2D' 1 



(52) 



1 _ Pj-+j 

m=l J J 

where Eq. ffl2|) was used, the subscript "1" denotes N — 1 and the subscript "gT stands for 
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"discrete." The second moment, or the mean-square FPT (MSFPT), is 



(Ti )d = t>J2 m \ m = T f^i =m^ + ^ (53) 



m=l 



On the other hand, the corresponding continuum results are 21] 



W. = (54) 

(n% - ^ (55) 

Note that while the results for the MFPT are always the same {(Ti)d — (Ti) c ), those for 
the MSFPT only coincide for Pj->j = 2/3, which corresponds to the optimal algorithm. 
Note also that the MSFPT for the LMC [Eq. ([53}] is the lowest for p^ = (the "trivial" 
ordinary algorithm). In fact, in this case the MSFPT is simply the square of the MFPT, 
i.e., the variance of the FPT is zero — the FPT is deterministic and is always equal to the 
time step of the algorithm. On the other hand, the optimal algorithm produces the correct 
variance of the FPT. 

In fact, the matching of the first and second moments of the FPT for N = 1 guarantees 
such matching for any other N and even for "asymmetric" FPT problems where the two 
walls are at different distances from the initial position. This can be seen from the following 
consideration. Consider a random walk in continuum space and map it onto a lattice walk 
by introducing sites at points x = ja and making a jump of the lattice walk from site j to 
site j ± 1 when the continuum walk visits site j ± 1 for the first time after visiting site j. 
The probability of any particular sequence of visited sites in such a walk will be the same as 
in the LMC algorithms (either the ordinary or the optimal one), since the only requirement 
for the correct sequence is that the walk be unbiased. This means that the probability 
P m to reach in m jumps a particular set of sites of interest starting from another specific 
site will always be the same for the walk obtained by direct mapping from the continuum 
walk and for any of the LMC approximations (where for the optimal LMC only the actual 
jumps are counted and not the steps where the particle stays put). On the other hand, the 
mean waiting times for a single jump are given by Eq. (152]) for the LMC algorithms and by 
Eq. ([MP for the walk obtained by mapping from the continuum walk; likewise, the mean- 
square waiting times are given by Eqs. (l53~l) and ( l55l) . respectively. Therefore, if it is known 
that the first-passage number of jumps is m, then the MFPT is m(Ti)d )C and the MSFPT 
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is m{m — l){Ti) 2 dc + m(T 2 )d tC , where the appropriate subscript d or c and the appropriate 
Pj^j are chosen depending on the walk. After averaging over all possible m we get for the 
MFPT 

(T) = Y,PnMTi) d ,c (56) 

and for the MSFPT, 

<T 2 > = £P m (m(m - lXH}^ + m^V). (57) 
m 

Given that P m are the same for all three walks and since (Ti)d — (Pi)c for all Pj->j, the 
MFPT for both LMC algorithms is always the same as for the walk obtained by mapping 
the continuum walk and thus the same as for the continuum walk. This is also true for 
the MSFPT for the optimal LMC algorithm, as in that case (T 2 )^ = (T 2 ) c , but not for the 
ordinary LMC, as in that case this equality does not hold. Instead, for the ordinary LMC, 
using Eq. (J57J) and Eqs. (|52J)-(|53j) with = 0, 

(T 2 ) = ^P m [m(m- l)(a 4 /4P 2 ) + m(a 4 /4P 2 )] = ^ £ m 2 P m = ^(™ 2 >, (58) 

m m 

where (m 2 ) is the mean-square first-passage number of jumps of the unbiased random walk. 
This is just the square of the (deterministic) interval between jumps (equal to r) times (m 2 ). 
For the FPT problem where the particle starts at site and the walls are at sites ±iV, (m 2 ) 



is 



22 1 



. 2v 5iV 4 -2iV 2 . 
(m 2 ) = , (59) 



so 

,4 



(^ 2 ) = I ^(5iV 4 -2iV 2 ). (60) 
If a = b/N, then 

(T 2 ) = T ^(5-2/iV 2 ), (61) 

which coincides with the continuum result 56 4 /(12P 2 ) in the limit N — > oo, but deviates 
for finite N. This is illustrated in Fig. |3]^a), where the comparison to the optimal algorithm 
is made; for the latter, as mentioned, the exact result is recovered for all N. Note that the 
data in Fig. [3] are obtained for b = 1, which gives a = 1/N, and so iV = 1/a is simply the 
inverse mesh step. 

In Fig. ((Sj)(b), we show the third moment of the FPT. These results have been obtained 
numerically, by iterating the corresponding master equations, which, of course, is much 
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FIG. 3: For the ordinary [Eqs. (PJ-©] and optimal [Eqs. f[26|> — (p8]) ] LMC algorithms, the first 
two moments [MFPT and MSFPT; panel (a)] and the third moment [panel (b)] of the FPT that 
would be observed in the simulations of ID diffusion with D = 1/2 starting at x = at t = and 
with the walls located at x = ±1 (i.e., b = 1), as a function of the inverse mesh step N. The data 
points are obtained by numerically iterating the corresponding master equations. For the optimal 
algorithm, the MFPT and MSFPT are independent of N and coincide with the continuum values 
(1 and 5/3, respectively, shown by the dotted lines). For the ordinary algorithm, the MFPT is 
likewise A^-independent, but the MSFPT data lie on the curve given by Eq. (161 f) (dashed line). For 
the third moment, the dashed lines are perfect fits to the data, and the dotted line represents the 
continuum value. 
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more accurate than would be achievable with actual LMC simulations. As expected, the 
results are now iV-dependent in both cases, but the deviation from the asymptotic value is 
much smaller for the optimal algorithm and, in fact, is barely visible already for N = 1. 
While we have not derived the analytical results (this would be very tedious to do), the 
equations given in the figure provide essentially perfect fits to the data, and we believe 
that these are exact (with any deviations due entirely to numerical errors). It is seen that 
while the deviation from the continuum value is 0(1/N 2 ) for the ordinary algorithm, it is 
0(1/N 4 ) for the optimal one, in full analogy with the higher moments of the spatial particle 
distribution (Sec. Ill B[) . where the leading term in the deviation from the continuum value 
likewise vanishes when the algorithm is optimized. 

We can also obtain numerically, likewise by iterating the corresponding master equations, 
the whole FPT distributions for different algorithms and compare them to the continuum 
result. First of all, the immediate shortcoming of the ordinary algorithm is that for even N, 
the particle can only reach the boundaries at even steps, while for odd N, it can only reach 
the boundaries at odd steps. This means that the first-passage probability will vary wildly 
alternating between zero and a nonzero value even for arbitrarily large N. To eliminate 
this, we plot only the nonzero values of the probabilities dividing them by two, in effect, 
averaging between the zero and nonzero values. 

The continuum result for the probability density r(t) of reaching the boundary at time t 



(or the first-passage rate) is (see, e.g., Ref. [23[) 



n . , 

r(t) = ^(-l) m (2m + 1) exp(-(2m + 1) \ 2 Dt/4b 2 ). (62) 

m=0 

This is a function that is exponentially small for small t, reaches a maximum and decays 
exponentially for large t (the thick solid line in Fig. HJ). 

Take for simplicity b — 1 and D = 1/2, which gives the MFPT (T) = 1. In this case, 
a = 1/N and the time step is 1/N 2 for the ordinary algorithm and l/(3iV 2 ) for the optimal 
algorithm. The comparison of the ordinary and optimal algorithms (with different degrees of 
discretization N) with the exact result is made in Fig.|H For N = 1, the FPT is deterministic 
and always equals 1 for the ordinary algorithm, as the particle always moves to site +1 or —1 
during the first step. For the optimal algorithm, on the other hand, the particle can move 
to one of the sites ±1 at any step with probability 1/3; therefore, the FPT distribution 
decays exponentially. While the maximum present in the continuum distribution is not 
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Time t 



FIG. 4: The continuum first-passage rate for ID diffusion with D = 1/2 starting at x = at 
t = and with the walls located at x = ±6 = ±1 (thick solid line), compared to the corresponding 
quantity for the LMC algorithms (symbols, in some cases connected by thin lines for guidance 
to the eye). This first-passage rate analog is calculated as the probability of first passage at a 
particular step divided by the time step of the algorithm. For the ordinary algorithm, the even 
steps for odd N and the odd steps for even N, at which the passage probability is always zero, are 
discarded and the rates for remaining time steps are divided by two. The inset shows the difference 
between the LMC rate and the continuum rate for the optimal algorithm with N = 2 (thick line) 
and for the ordinary algorithm with N = 6, 8, 10 (thin line), as indicated by the labels next to the 
curves. 

reproduced, the continuum distribution is matched quite closely for large t. Already for 
N = 2 the optimal algorithm reproduces the continuum curve very accurately. With the 
ordinary algorithm, a similar accuracy is not achieved until N rj 8 — 10, as the inset of 
Fig. H] shows, where the difference between the LMC results and the continuum values is 
plotted. Very good results at large t obtained with the optimal algorithm would not be 
surprising, as it is designed specifically to be as accurate as possible in this case; however, 
we see that it also works very well for moderate t. In terms of the computational effort, the 
ordinary algorithm with N = 8 takes N 2 = 64 steps on average to reach the walls, whereas 
the optimal algorithm with N = 2 takes 3iV 2 = 12 steps on average; thus the speedup is a 
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factor of at least 5. The advantage is even larger when the "numerically exact" approach 



based on solving numerically the master equations [20[ is used, since decreasing the mesh 
step increases the number of such equations and not only the number of time steps. 

As mentioned above, the accuracy of the optimal algorithm is expected to be particularly 
good at large times. In the limit of large t, the first-passage rate decays exponentially: 

r(t) ~ 7exp(-/3i). (63) 

It is therefore convenient to compare the algorithms by finding the rate of the exponential 
decay (3 and the prefactor 7 and comparing to the continuum values. By looking at the first 
term of Eq. fl62|) . for D = 1/2 and b = 1 we get the continuum decay rate of /3 C = 7r 2 /8 and 
the prefactor of 7 C = tt/2. For the LMC algorithms, these quantities can be derived as well 
(see Appendix). For the ordinary algorithm, the decay rate is 

/3 ord = -iV 2 lncos(^), (64) 



the prefactor is 



7 ord = iVtan(^). (65) 



For the optimal algorithm, they are respectively 



/3 opt = -3A^ln( jj + ^cos(^) ) (66) 



and 



^sin(7r/2iV) 

7opt 2/3 + (l/3)cos(7r/2A0' 1 ' 

These quantities are plotted in Fig. [5j For the optimal algorithm, even for N = 1 the results 
are quite close to the continuum value (/3 opt = 1.216 ... vs. tt 2 /8 = 1.233 . . . for the rate 
and 7o P t = 1.5 vs. tt/2 = 1.570. . . for the prefactor). For the ordinary algorithm, similar 
accuracy is achieved for N = 4 — 5; note that for N = 1 the decay rate diverges, as the FPT 
is deterministic. For N = 2 the results for the optimal algorithm are nearly indistinguishable 
from the continuum values on the scale of the plot and are closer than those for the ordinary 
algorithm with N = 10. Note that the long-time approximation is applicable even for quite 
small t: for example, for t = 1 (equal to the MFPT), the full sum in Eq. (1621) is 0.457365 . . . 
and its first term is 0.457436 . . ., so the difference is only f» —7 x 10~ 5 . 
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Inverse mesh step N 

FIG. 5: The large-time rate of exponential decay [panel (a)] and the corresponding prefactor [panel 
(b)] for the first-passage rate to the walls located at x = ±1 starting at x = in simulations of ID 
diffusion using the ordinary (open circles) and the optimal (filled circles) LMC algorithms. The 
curves are the analytical expressions [Eqs. (|64p - (|67p ]. The continuum limit values are indicated by 
the dotted lines. 



III. TWO DIMENSIONS 



In 2D, the continuum unbiased diffusion equation is 



dn(x,y,t) _ D / d 2 n(x,y,t) | d 2 n(x,y,t) ' ^ 



dt \ dx 2 dy 2 
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Its general solution can be written in a form similar to Eq. (|T5|) : 

/OO POO 
/ C(k x , k y ) exp(ik x x + ik y y - a c (k x , k v )i) dk x dk y . (69) 
-oo J — oo 

The dispersion relation is analogous as well: 

a c (k x ,k y )=D(k 2 x + k 2 y ). (70) 

As for the design of the corresponding LMC algorithm, we first note that different types 
of lattices can be chosen. We consider the simplest choice, the square lattice, with the lattice 
constant still denoted a. The lattice sites can be denoted by pairs of integer numbers, 
If only moves along the Cartesian axes are allowed, there are five different probabilities of 
moves: p^j^^ + ij) = p +x and P(i,j)^(i-i,j) = P-x for the moves along the x axis in the 
positive and negative directions, respectively; p^-^jj+i) = p +y and p^j^^j-i) = p- y for 
the moves along the y axis; and the probability of staying on the same site, P(i,j)-^(ij) = Po- 
The master equation for the mean particle number nun at a particular site (J, I) is 

n m {t + T ) = Pa n U,i) (*) + P+xTiQ-1,1) (t) + p- x n(j +1 j) (t) + p+ynuj-!) (t) + p- y n(j )l+1 ) (t) . (71) 

The general solution of the system of master equations for all sites is (again, similarly to 
ID) 

/TV /a pn/a 
/ C(k x , k y ) exp(ik x aj + ik y al - a d (k x , k y )t)dk (72) 
•7r/o J —n/a 

with 

ad{k x ,k y ) = -- \n(p + p +x exp(-ik x a) + p_ x exp(ik x a) + p +y exp(-ik y a) + p_ y exp(ik y a)) . 

(73) 

As in ID, the goal is to choose the probabilities and the time step r so that a c and ad match 
as closely as possible for small k x and k y . While it is possible to match all coefficients of the 
Taylor expansion up to and including the quadratic terms (which, as in ID, guarantees the 
correct diffusion rate), matching all quartic terms, like we did in ID, is obviously impossible. 
This is because the form of Eq. (173")) is such that there are always terms oc k 2 x and cx ky, 
but no quartic term oc k x k 2 under the logarithm. When the logarithm is expanded, this will 
necessarily produce such a quartic term. At the same time, no quartic terms are present 
in the expression for a c [Eq. (ITU]) ] that we need to match. This reasoning is valid for any 
set of moves, even if moves to second neighbors and more distant sites are present, as 
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long as all moves are along one of the axes. On the other hand, moves along both axes 
simultaneously will produce a term cx k 2 ky under the logarithm with a coefficient that can 
be chosen so that in the final expansion this term vanishes. The simplest moves of this 
type are the four "diagonal" moves whereby the particle moves by one lattice constant along 
the x axis and simultaneously along the y axis (Fig. E]). This introduces four additional 
probabilities: p^^+ij+i) = p +x , +y , P(i,j)^(i+i,j-i) = P+x-y, P(i,j)-Ki-i,i+i) = P-x,+y, an d 
= p- x -y The new master equation is 

n(j ; i)(t + r) =p n {j:l) (t) + p +x n {j „ ltl) (t) + p- x n {j+1)l) {t) + p +y n (t) + p- y n {j j +1) {t) 

+ P+x,+ y n(j-i t i-i){t) + p +x - y n (j _ lt i +1) (t) 

+ p_ :E)+2/ ri( i+1)i _i)(t) +p_ (B ,_yn( J - + i,i + i)(t). (74) 
The dispersion relation now becomes 

a d (k x ,k y ) = — ln(p + P+x exp(-ik x a) + p_ x exp(ik x a) + p +y exp(-ik y a) + p^ y exp(ik y a) 

+ p+ x ,+ y exp(—ik x a — ik y a) + p+ x - y exp(—ik x a + ik y a) 
+ P- x ,+ y exp(ik x a — ik y a) + p~ x -y exp(ik x a + ik y a)). (75) 

Matching the k° term, 

p + P +x + P- x + P +y + P-y + P+x,+y + P+x-y + P-cc,+y + P-x-y = 1, (76) 

which is again just the normalization statement. Matching the linear terms, 

P+x + P+x,+y + P+x-y ~~ P-x ~ P-x,+y ~ P-x-y = 0, (77) 
P+y + P+x,+y + P-x,+y ~ P-y ~ P+x-y ~~ P-x-y = 0. (78) 

These ensure that the average velocities are zero in both directions (no bias). Taking 
Eqs. (!76|) - (l78l into account, we can expand under the logarithm to obtain 

a d (k x , k y ) = — - ln(l - (C xx /2)a 2 k 2 x - (C yy /2)a 2 k 2 y + C xy a 2 k x k y - i(C xxy /2)a 3 kX 

- i(C xyy /2)a 3 k x k 2 y + (C xx /2A)a A k x + {C yy j2£)a% - (C xy /6)a 3 k 3 x k y 

- (C xy /6)a 3 k x k 3 y + (C xxyy /A)a A k 2 x k 2 + 0(k 5 )), (79) 
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where 

C xx — P+x + P-x + P+x,+y + P+x-y + P-x,+y + P-x-y = 

2(j) +x + P+ x ,+y + P+x-y), (80) 

Cyy = P+y + P-y + + P+x-y + P-rr,+y + P-x-y = 2{p+ y + + P-x,+y), (81) 

Cxj/ = — + P+x-y + P-x,+y ~~ P-x-yi (82) 

C X xy "I - P+x—y P—x,+y P—x,—yi (83) 

C'xyy P+x,+y P+x,—y P—x,+y P—x,—yi (84) 

C X xyy = P+x,+y + P+x,-y + P-x,+y + P-x-y (85) 

Expanding the logarithm, we now get 

&d(k x , ky) = — ((C xx /2)a k x + (C yy /2)a k y C xy o, k x k y + i(C X xy/2)Q?k x k y 

•'(r,, w; 2)r,^-,i; + (C«,/8 - C M /24)o^ + (C 2 J8 - C yy /2A)a% 
~^(C xy /6 — C xx C xy /2)a 3 k x k y + (C x . y /6 — C yy C xy /2)a?k x k y 
HCxxCyy/A - Cxxyy/^ + C 2 xy /2)a A k 2 x k 2 y + 0(k 5 )). (86) 

Comparing to Eq. (ITOl) . we obtain the following equations: 



C X x°? /2t 


= D, 


(87) 


CyyQj /2T 


= D, 


(88) 


C xy /T 


= o, 


(89) 


C X xy /T 


= o, 


(90) 


C x yy/T 


= o, 


(91) 


(C 2 xx /8 - C M /24)/r 


= o, 


(92) 


(C*,/8 - C w /24)/t 


= o, 


(93) 


(C xy /6 — C xx C xy /2) /t 


= o, 


(94) 


(C X y/6 — C yy C X y/2)/T 


= o, 


(95) 


-C XX yy/A + Cl y /2)/T 


= 0. 


(96) 



From Eqs. (152)1 and (Pjh C xx = C yy = \ (C xx = and C yy = are not solutions, as in that 
case r = 0). From this, using Eq. ( 1H71) or ( |88l) . we can compute the time step: 

' = 55- (97) 
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Interestingly, this optimal time step is equal to the one we obtained previously for the ID 
problem [Eq. fl28l) ]. From Eq. f )89|) . C xy = 0. By inserting the values of C xx , C yy and C xy in 
Eq. (jnS]), we obtain C xxyy = 1/9. From Eqs. fl9U]) and (J9JQ), C xxy = C xyy = 0. Equations (ED 
and f l93|) are then satisfied automatically. Thus all coefficients C are known. Using the 
expressions for C [Eqs. (j80p — (I85p ]. we obtain linear equations for the probabilities. Together 
with Eqs. (!77|) and ( !7H|) . this forms a linear system of eight equations with eight unknowns 
that has a unique solution, 

P +X = P- X = P+y = P-y = 1/9, (98) 

P +Xj+y = P+x-y = P-x,+y = P-x-y = 1/36. (99) 

Finally, from Eq. (ITB"j) . 

Po = 4/9. (100) 

The set of parameters given by Eqs. (I97j) - (ll00p defines an unbiased 2D LMC algorithm 
that is "optimal" for the particular set of moves considered here, since it is the only set 
of parameters that reproduces the continuum dispersion relation up to 0(k 5 ), the same 
accuracy as for the optimal ID algorithm considered in Sec. [TT1 [In both cases, the k 5 terms 
were not required to vanish explicitly, but are zero automatically, as are all other odd-order 
terms.] The set of moves is itself optimal in the sense that it is the smallest set of shortest- 
ranged moves that ensures such accuracy. We note that the probability of staying put is 
reduced from 2/3 to 4/9, and that the diagonal moves have small but finite probabilities. 

Note that the probabilities (I9"8"|) - (ll00p can be obtained as products of ID optimal probabil- 
ities [Eqs. ([52]) and (I27p ]: the probability of every 2D move is the product of the probabilities 
of the ID moves that it "consists of." For example, the 2D diagonal move (-He, +y) can 
be thought of as consisting of two ID moves, j — 1 — > j, one in the x direction and one 
in the y direction, and indeed, p+ x ,+ y = Pj-i^j x Pj-i^j- Likewise, a move along one axis, 
say, —x, in 2D LMC can be represented as a move by one site in that direction and staying 
put in the orthogonal direction, which gives the correct equality p_ x = Pj+i^j X Pj^j. For 
the probability of staying put, by similar reasoning p Q = Pj^j x Pj^j, which is also correct. 
Thus the 2D optimal algorithm can be considered as the direct product of the ID optimal 
algorithm with itself, in the sense that if the 2D probabilities are arranged in a matrix, this 
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(i-lj-l) O'.J-l) 



FIG. 6: A schematic of the general 2D LMC algorithm considered in this paper. Only moves to and 
from site (i,j) are shown. In addition to the traditional axial moves (arrows of medium thickness), 
diagonal moves (thin arrows) are introduced; also, there is a probability to stay put (represented 
by the thickest arrow). Notation for the probabilities of the moves is given. 



matrix is the direct product of the vectors of the ID probabilities: 



P—x,—y P—x P—x,+y 
P-y PO P+y 
\P+x-y P+x P+x,+y J 



Pj+l-H X Pj • I \j Pj ■ I V X Pj './ Pj+l->j X Pj • *j 

Pj •./ x Pj+l->j Pj './ x Pj •./ Pj :i x Pj 



Pj+l-*j I'i •; Pj- 



(101) 



Pj+1-,3 
Pj^j 

The fact that the direct product algorithm should work as well in 2D as its "fac- 
tors" do in ID can be seen directly from Eq. (1751) : if the probabilities entering that 
equation are the products of the ID probabilities, then the expression under the log- 
arithm can be represented as a product of two expressions, f{k x )f(k y ), where f(k) = 
Pj^j + Pj-i^j exp(— ika) + Pj+i->j exp(ika) is the expression under the logarithm in the 
ID analog of Eq. fl75|) . Eq. f|T9|) . But if the probabilities in f(k) are chosen optimally, 
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then, according to Eq. (|T9|) and the considerations following it, ln(/(fc)) = —rDk 2 + 0(k 6 ), 
therefore in Eq. ( 175]) 

a d (k x , k y ) = -i ln[/(*,)/fo)] = [Dk£ + 0(*£)] + [DA* + 0(k% (102) 

which matches the continuum dispersion relation [Eq. (ITU] up to and including the 4th order 
terms. (Note, by the way, that there are no cross terms containing both k x and k y in any 
order and also, as in ID, no odd-order terms.) Note that this consideration assumes that 
the time step in 2D is the same as in ID, which is indeed the case. 

Moreover, the statement that the direct product algorithm is as good in 2D as its factors 
are in ID is true whenever the solution of the original 2D equation can be represented as 

n(x, y,t) = J C(k x , k y )g kx (x, t)g ky (y, t)dk x dk y , (103) 

where g k (x,t) are some solutions (depending on the parameter k) of the ID equation for 
which the ID algorithm works. This can be checked by substituting gk x (%, t)gk y (y, t) into the 
master equation ( 1741 [much like we did in Eq. (132]) ]. which makes both sides of that equation 
products of the corresponding sides of the ID master equation (TTTjl . in which g^ x and g ky 
are likewise substituted. In this sense, we could have immediately "guessed" the correct 
2D algorithm by looking at the solution of the 2D diffusion equation [Eq. ( 169]) . from which 
gk(x,t) = exp(ikx — Dk 2 t)}, without the lengthy derivation that followed. That derivation 
was still useful, however, as it proved both the uniqueness of the solution given the set of 
moves and the fact that the diagonal moves were necessary. 

Finally, note that from the fact that the 2D algorithm is the direct product of the ID 
optimal algorithm with itself, it follows that projecting all moves of the 2D optimal algorithm 
along the x axis and keeping the time step the same produces the ID optimal algorithm. 

IV. THREE AND MORE DIMENSIONS 

In 3D, first of all, we can use the same arguments as in 2D to prove that the 3D algorithm 
constructed from the "direct product" of the ID algorithms will reproduce the dispersion 
relation with the same precision, i.e., up to 0(k 5 ). This 3D algorithm involves four types of 
moves: in addition to moves along one of the Cartesian directions, diagonal moves along two 
directions and staying put, the new type of move in 3D is a move along all three directions 
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simultaneously. The probabilities of the moves are 

po = (2/3) 3 = 8/27, (104) 

p ±x = p ±y = p ±z = (2/3) 2 x (1/6) = 2/27, (105) 

P±x,±y = P±y,±z = P±x,±z = (2/3) x (1/6) 2 = 1/54, (106) 

P±*,± y ,±z = (1/6) 3 = 1/216, (107) 
and the time step is the same as in ID and 2D: 



- ^_ 
' ~ 6D' 



(108) 



However, unlike in 2D, the direct product algorithm is not the only way to satisfy the 
dispersion relation up to 0(k 5 ). For instance, there is no need to introduce moves simul- 
taneously along all three directions: just like the moves along two directions are used to 
adjust the term oc & 2 /c 2 , the moves along three directions would be useful to adjust the 
term oc k x kyk^, but this term is 0(k 6 ), and the terms of this order cannot all be made to 
vanish simultaneously anyway, since there are not enough different moves (i.e., not enough 
free parameters). Therefore, we can retain just the moves that we used to optimize the 2D 
algorithm, i.e. those along either one or two directions, and the resulting equations are very 
similar. In particular, the analogs of Eqs. ([75]) . ([75]) and (JS5|) for ct^ are obtained simply 
by adding the analogous terms involving the z axis. This is also true for the analogs of 
Eqs. (ITS"]) . (ITT]) and (ITS]) , to which a similar equation obtained from equating the coefficient 
of k z to zero is added. Expressions for the coefficients C in terms of the probabilities p [the 
analogs of Eqs. (!80]) - (!85l) ] as well as those for the new coefficients involving the z axis (C zz , 
C xz , C yz , C xxz , C yyz , C xzz , C yzz , C xxzz , Cy yzz ) are built by analogy as well. The equations 
for C [Eqs. (18T]) - (19"6"]) ] are still valid, and the analogous equations involving the z axis are 
added to them. From these equations, C xx = C yy = C zz = 1/3, C xxyy = C xxzz = C yyzz = 1/9 
and all other C are zero; the time step is again the same as in ID and 2D, as well as in the 
direct product algorithm in 3D, i.e., it is given by Eq. (I108p . Finally, for the probabilities 
we obtain: 

P +X = P- X = P + y = P^y = P +z = P- Z = 1/18, (109) 

P+x,+y = P+x-y = P-x,+y = P-x-y = P+x,+z = P+x-z = 

P-x,+z = P-x-z = P+y,+z = P+y,-z = P-y,+z = P-y-z = 1/36, (110) 

Po = 1/3. (HI) 
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Projecting the moves of either this algorithm or the direct product algorithm [Eqs. ( I104p - 
(jl08p ] onto the x axis produces the ID optimal unbiased algorithm and projecting onto the 
xy plane produces the 2D optimal unbiased algorithm. 

We now have two different "optimal" LMC algorithms for 3D. The second algorithm is 
the only possible one without moves in three directions; however, if these moves are allowed, 
the first algorithm is but one example. This "direct product" algorithm, in a way, has an 
advantage that some of the sixth- and even higher-order terms in the dispersion relation 
vanish; in fact, it follows from a consideration analogous to that of the previous section that 
all cross terms containing at least two of the three variables k x , k y and k z are zero in all 
orders. However, since other sixth-order terms (k®, ky and k\ ) do not vanish and given that, 
generally speaking, all sixth-order terms are expected to be of about the same magnitude, 
there is no guarantee that the "direct product" algorithm is more accurate than either the 
algorithm without moves in three directions or other possible algorithms where such moves 
are allowed. This may be problem-dependent. Note that fifth-order terms (as well as all 
other odd-order terms, in fact, even those odd in any one of the variables) vanish in both 
algorithms. On the other hand, the advantage of the algorithm without moves in three 
directions is that it is somewhat more short-ranged and slightly easier to program, as there 
are fewer moves. 

We note that in the 3D case, the time step of either of the two "optimal" algorithms we 
have considered, r = a 2 /6D, is also the time step of the non-optimal but straightforward 
"ordinary" algorithm where the particle moves at each step along one of the six Cartesian 
directions. This is unlike in ID and 2D, where the analogous straightforward algorithms 
have a larger time step and thus are more computationally efficient (as we discussed in 
detail for ID). Note, however, that compared to the straightforward algorithm, we now use 
an extended set of moves including longer-range diagonal moves that allow an even larger 
time step when optimal accuracy is not required. That a larger time step is possible with 
the expanded set of moves is particularly obvious in 2D, where the algorithm with diagonal 
moves only is equivalent to the straightforward algorithm with the mesh step of v2 times 
the original mesh step (and the axes rotated by 45°) and thus the corresponding time step 
is twice that of the original straightforward algorithm. 

Extending these results to a space of arbitrary dimensionality using a hypercubic lattice 
is straightforward. First of all, "direct product" algorithms can be constructed, with the 
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same properties as in 2D and 3D. On the other hand, the argument why simultaneous moves 
in more than two directions are unnecessary still applies (however, avoiding such moves may 
be impossible for other reasons, as we explain below). To construct an algorithm with moves 
in one direction (axial moves) and two directions (diagonal moves) only, equations analogous 
to Eqs. f )87p - (l96l) can be written for any dimensionality, and the solutions are always the 
same: all coefficients C with two equal indices (e.g., C xx ) are always equal to 1/3, all such 
coefficients with two pairs of equal indices with no equality between pairs (e.g., C xxyy ) are 
equal to 1/9, and all other C are zero. In d dimensions, 

C xx = 2 Pl + A(d - l)p 2 , (112) 

C xxyy = 4^2, (H3) 

where p\ is the probability of any particular axial move and p 2 is the probability of any 
particular diagonal move. This gives 

P. = ^, (114) 

n = (H5) 

It is easy to check that this works for d = 1, 2 and 3 (except for d = 1, there are no diagonal 
moves and so p 2 is meaningless). For d > 3, though, p\ is only non-negative for d = 4. In 
this case, p\ = 0, so the optimal algorithm contains no axial moves, only diagonal ones. 
Since there are ( 2 ) = 6 pairs of axes and thus 2 x 2 x 6 = 24 different diagonal moves, 
the total probability of a diagonal move is 24/36 = 2/3, so the probability to stay put is 
Po = 1/3, as in 3D. The time step is still a 2 /6D. Note that this is now larger than in the 
straightforward algorithm with only axial moves and zero probability to stay put; but again, 
a larger time step is only possible because we use longer-range diagonal moves. For d > 4, 
it is impossible to make all k 4 terms in the dispersion relation zero with the set of moves 
considered. This may be possible using even longer-ranged moves (e.g., moves to second 
neighbors along axes) and/or simultaneous moves in three or more directions. In particular, 
the "direct product" algorithm is always an option in space of any dimensionality. 

V. TREATMENT OF IMPENETRABLE (REFLECTING) BOUNDARIES 



So far, we have treated diffusion in a homogeneous medium without any boundaries. Of 
course, this is a trivial situation and the really interesting cases are those where some spatial 

30 



inhomogeneities (e.g., in the form of obstacles) are present. In this section, we examine 
how to take impenetrable "reflecting" boundaries into account in LMC algorithms. Such 
boundaries are introduced into the continuum diffusion problems via boundary conditions, 
which in the unbiased case reduce to the expression 

b-^n = 0, (116) 

where b is the unit vector normal to the boundary. Since the particle flux is proportional to 
the concentration gradient \^n, Eq. (jllfip ensures that there is no particle flow across the 
boundary. 

A. One dimension 

We consider the case of a single boundary (or wall), which means that particles cannot 
go beyond a certain point x = Xj, on the line; without the loss of generality, we choose 
Xb = and assume that the region x > is accessible (Fig. [7]). Let the site nearest to the 
wall have index (the filled square in Fig. [7]). Since this site has a neighbor only on the 
right, moves from the left to that site are impossible, and the master equation for the mean 
particle number at that site has to be modified: 

n (t + r) =po_H)»io(t)+Pi-+oni(t). (117) 

All other sites (which we refer to as bulk sites) have two neighbors and we assume that 
the corresponding master equations, including all the probabilities that are involved and 
the time step, remain as derived previously for free space (Fig. [7J. Given that assumption, 
the probabilities entering Eq. (JUTJ, Po->o an d Pi-»o, can be determined uniquely from the 
condition that the probabilities of all moves leaving a given site (including the move to the 
same site, j — > j) should sum up to 1 (the normalization condition). In particular, for site 
we have 

Po->o + Po->i = 1. (118) 

However, po->i enters the master equation for site 1, which, by our assumption, remains the 
same as in free space, so Pcm-i is the same as Pj-i-^j in free space, which determines 

Po^o = 1 - Pj-i-*j = Pj+i^j + Pj-*j- (119) 
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FIG. 7: Treatment of a reflecting boundary in the ID optimal LMC algorithm. The inaccessible 
region is shaded. The site closest to the boundary is the boundary site (filled square); all other 
sites are bulk sites (filled circles). The open square is a fictitious site in the inaccessible region 
useful in the analysis of the algorithm, as described in the text. The optimal placement of the 
boundary site is at the distance d = a/2 from the wall. Only the probabilities of the two moves 
leading to the boundary site (thick arrows), including the move from the boundary site to itself, are 
allowed to change (although pi^o does not in the end); the probabilities of all other moves (thin 
arrows) are assumed to be the same as in the bulk from the outset. Notation for the probabilities 
of the moves used in the text is given. The move to the left from the boundary site (dashed arrow) 
is rejected and contributes to the probability of staying put. 

Likewise, for site 1 



Here is involved in the equation for site 1 (and thus equals the "bulk" Pj->j) and pi_>2 is 
involved in the equation for site 2 (thus equals the "bulk" and so pi^o is determined 

uniquely as 



The time step should necessarily be the same as for bulk sites. This is because all moves out 
of a site must have the same time step associated with them, and for site 1, while a move 
1 — > leads to a boundary site and thus is considered a "boundary" move, a move 1 — > 2 is 
a "bulk" move. So the time step for the "boundary" moves is uniquely determined. 

Then, for instance, for the algorithm for unbiased diffusion optimized for accuracy, 



= 1. 



(120) 



Pi^o = 1 - Pj^j ~ Pj-l^j = Pj+l^j- 
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Eqs. (]26D (J28D lead to 



Po^o = 5/6, (122) 

Pi^o = 1/6, (123) 

Po^i = 1/6, (124) 

r = a 2 /6D. (125) 

The probability po->i is a bulk probability and is given here for completeness. Note that the 
probability of moving to the left from site 1 is the same as from any site with index j + 1 > 1 
(Pi ->o — Pj+i->j)- However, the probability of staying put increases by 1/6. This quantity 
1/6 can be interpreted as the probability of the move to the left from site that is now 
forbidden [in fact, this interpretation follows directly from Eq. ( 11191) ]. In other words, the 
algorithm is: for the particle in any site, including the boundary site, select the next move 
according to the bulk probabilities, but if the boundary is crossed, reject the move and stay 
on the same site. 

We can check explicitly that using Eqs. H122|) — fll25[) in the master equation fl 1 1 7[) indeed 
provides the best approximation of the continuum equation. In ID, the boundary condition 
[Eq. 0116p ] becomes 

f)n\T I i 

0. (126) 



dn(x, t) 



9x x=0 

The general solution of the diffusion equation [Eq. (jSJ)] with this boundary condition can be 
written as 

;>oo 

n{x,t)= C(k)cos(kx)exp(-Dk 2 t)dk, (127) 
Jo 

which is just the general solution in free space [Eq. f)15p ] with the restriction C(—k) = C(k) 
arising from the reflecting boundary condition f)126p [which, together with the original con- 
dition C(—k) = C*(k), means that C(k) is real]. Note that because of broken translational 
symmetry, generally speaking, the solution of the discrete master equation can no longer be 
written in the same form we used in the bulk [Eq. ffT8l) ]. Therefore the dispersion relation 
cannot be introduced and we have to use the alternative method based on inserting the 
integrand of the solution of the continuum equation [Eq. (11271) ] into the master equation 
and requiring that it turns into an identity as accurately as possible for small k [see Eq. 
and the accompanying discussion]. 
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Even though the master equations for the bulk sites involve only the bulk probabilities 
that we have already determined in Sec. [Til h is still instructive to sketch the consideration 
for these sites using the new form of the solution [Eq. (11271) ]. as there are some subtleties 
compared to Eq. (I31Z1) . Inserting the integrand of Eq. (I127p into Eq. ffT7|) and dividing both 
sides by C{k) exp(-DkH), we get 

? 

cos(kxj) exp(— Dk 2 r) ~ Vj-^j cos(kxj) +p 3 _ 1 _ > _ 7 - cos(k(xj — a)) cos(k(xj + a)), (128) 

where Xj is the position of site j. Note, however, that expanding straightforwardly in a series 
in k would not be correct, since Xj can be arbitrarily large. We did not face this problem 
in Sec. HJJ, since dividing by exp(ikaj) eliminated all dependence on j. However, ka can still 
be assumed small, so the correct approach is expanding cos(k(xj + a)) in series in ka: 

cos(k(xj ± a)) = cos(A^)(l - k 2 a 2 /2 + fc 4 a 4 /24 + . . .) =f sm(kxj)(ka - k 3 a 3 /6 + ...). (129) 

The exponential on the left-hand side of Eq. (11281) can be expanded as normal. This gives 

cos(kxj)(l - Dk 2 r + D 2 k 4 r 2 /2 + . . .) ~ p^j cos(kxj) 

+ {pj-i^j + Pj+i-^j) cos(kxj)(l - k 2 a 2 /2 + k A a 4 /2A + . . .) 
+(pj_i_ >J - — pj + i^,j) sm(kxj)(ka — k 3 a 3 /6 + ...). (130) 

This should be satisfied for any j, thus for any Xj, which immediately gives Pj-i^j = Pj+i^j- 
Then after dividing both parts by cos(kxj) the dependence on Xj is eliminated, and we end 
up with a series in k on both sides. Matching the coefficients gives the equations for the 
parameters whose only solution is given by Eqs. f[26]) — (]28|) . as expected. 

To obtain the probabilities for the boundary site 0, we need to repeat the same procedure 
for that site. Note, however, that in addition to the probabilities, there is one other free 
parameter: the distance d from the boundary to site 0. In free space without boundaries the 
problem is translationally invariant, so displacement of all sites by the same amount does 
not influence the results, and we assumed, without any loss of generality, that site had 
coordinate Xq = 0. However, now that the translational symmetry of the problem is broken 
by the presence of the boundary, the coordinate Xq = d of site (Fig. [7]) is an important 
parameter. The coordinate of an arbitrary site j is thus Xj — d + aj. Then, inserting the 
integrand of Eq. (I127p in Eq. fl 1 1 T|) and dividing both sides by C{k) exp(—Dk 2 t), we get 

1 

cos(kd) exp(—Dk r) ~ Po^o cos(kd) + pi-^o cos(/c(d + a)). (131) 
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Dividing by cos(kd) and using r = a 2 /6D, as in the bulk, 

. 9>9 . . ? cos(k(d + a)) 

exp(-a 2 /c 2 /6) ~ p ^o + Pi^o ^ (fcrf) - ■ (132) 

Note that d ~ a, so fcrf is small, and the Taylor expansion can be done normally, without 
any complications. Expanding up to 0(/c 4 ), 

l-a 2 k 2 /6+a A k 4 /72 = (p ^ +p 1 ^ )-p 1 ^ {a 2 k 2 (l/2+d/a)-(a 4 k A /l2)(l/2+2d/a-4:(d/a) 3 )]. 

(133) 

This equality is only satisfied when the probabilities are given by Eqs. fll22l) - fll25p and 



d = a/2. (134) 

The meaning of the result for d and the algorithm as a whole can be understood as follows. 
On the one hand, any solution of the continuum diffusion problem with a boundary is also 
a solution of the problem in free space, without a boundary, if it is continued symmetrically 
beyond the boundary. On the other hand, there is a similar correspondence between the 
discretizations of these problems as well. Indeed, looking at the corresponding master equa- 
tions, the master equations for the bulk sites [Eqs. (IT?]) . (I2"6]) - (T2"8]) ] are, of course, identical 
for these problems, and the equation for the boundary site [Eqs. ( 1117p .( fT22l - (ll25l) ] in the 
problem with the boundary is also equivalent to the equation for a bulk site, if we introduce 
a fictitious site (numbered —1) to the left of the boundary (the open square in Fig. [7J) and 
assume that the mean particle number at that site is always equal to that at the boundary 
site = no). This matches the correspondence between the continuum problems, if the 
sites and —1 are located symmetrically with respect to the boundary, which immediately 
gives d = a/2. Note also that because of the equivalence between the master equations for 
the boundary and bulk sites, the solution of the system of all master equations can, in fact, 
be represented as 

n j(t) — / C(k) cos(kaj) exp (— ad(k)t) dk. (135) 



Jo 

[same as Eq. ffT8]) with C(k) real], even though, as mentioned, this is not possible in general 
for an arbitrary set of parameters p and d. Here j runs from to oo, i.e., it includes both 
the boundary site and the bulk sites. The dispersion relation ad(k) is, of course, exactly the 
same as in infinite space. 

Another situation of interest is diffusion in a finite interval bounded by two reflecting 
walls. Since in the consideration above for a single wall the effect of that wall on the 
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algorithm is local, only modifying the probabilities for a single boundary site, it should 
be possible to consider the two walls independently whenever there are at least two sites 
between them, with the modifications of the probabilities at each of the two boundary sites 
the same as in the case of a single wall. Introducing the fictitious sites behind the walls and 
repeating the above consideration leads to the same result. Of course, this assumes that the 
distance between the walls is a multiple of the mesh step a, so that it is possible to place 
both boundary sites at distance a/2 from the respective wall. This restricts the possible 
choices of a given the distance between the walls. 

B. Two dimensions 

In 2D, the situation is more complicated, especially since diagonal moves (some of which 
would cross the boundary) are involved. It is no longer possible to deduce the probabilities of 
the moves into the boundary sites based on the normalization conditions for the probabilities 
alone. We restrict our consideration to the case of a perfectly flat infinite boundary parallel 
to one of the Cartesian axes. By analogy with ID, the case of two parallel infinite boundaries, 
with the motion of the particle confined to the space between these boundaries, is easy, as 
each boundary can be considered independently. A more complicated and interesting case 
of finite and/or curved obstacles and boundaries is touched upon briefly in Sec. IVHI 

Consider specifically the case of a planar boundary located at x = so that the half-space 
x > is accessible (Fig. [8]). The boundary sites form a column at x = d (the filled squares 
in Fig. [8]) and have indices (0,1). For a site (J, I), the coordinates are Xj = d + aj, yi = al 
(there is still translational symmetry in the y direction, so we can choose arbitrarily y$ = 0). 
The master equations for the boundary sites are 

n( ,i)(t + r) =p' n m (t) + p' + „n (0ii _i)(t) + p'_ y ri(p tl+1) (t) + p'_ x n {U) (t) 

+ p'- x ,+ y n(i,i-i)(t) +p'- x - y n(i,i+i){t), (136) 

where the primed probabilities correspond to moves ending on boundary sites and may differ 
from the bulk (unprimed) probabilities, whose values remain as determined in Sec. III. 

The boundary condition H126[) is still true (except for the fact that the particle concen- 
tration n now depends on y as well), and the general solution of the diffusion equation with 
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FIG. 8: Treatment of a reflecting boundary in the 2D optimal LMC algorithm. The diagram only 
shows the moves to and from one of the boundary sites [(0, j)]. The inaccessible region is shaded. 
The sites closest to the boundary are the boundary sites (filled squares); all other sites are bulk 
sites (filled circles). Open squares are fictitious sites in the inaccessible region useful in the analysis 
of the algorithm, as described in the text. The optimal placement of the line of boundary sites 
is at the distance d = a/2 from the wall. Only the probabilities of the moves to the site (thick 
arrows) are allowed to change compared to the bulk sites. Notation for the probabilities of the 
moves used in the text is given. The move to the left from the boundary site (dotted arrow) is 
rejected and increases the probability of staying put, but the diagonal moves across the boundary 
(dashed arrows) are projected along the boundary increasing p'± compared to p± y . 

this boundary condition is 

/*oo /*oo 

n(x,y,t) — / / C(k x , k y ) cos(k x x) exp[ik y y — D(kl + ky)t]dk x dk y . (137) 

J k x =0 J k y =— 00 

Inserting the integrand in Eq. (11361) and remembering that Eq. (1971) for r should remain 
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valid, we obtain after dividing by C(k x , k y ) cos(k x d) exp[ik y al — D{k 2 x + ky)t\: 
exp(-a 2 (^ + ^)/6)~pi + p ' +y e~'"'" + p'_/ k '° 



+ (pU + + pU,-,^') C ° S ^(M) a>) ' (138) 

After expanding up to 0(/c 4 ) and equating the corresponding coefficients, we get: 

for A; : 

Pq + V+y + P~y + PLx + P- X ,+y + pL- -„ = lj (139) 



for A; 1 ,: 



for ky\ 



ia(P+y ~ P-y + P'-x,+y ~ P-x-y) = ^ ( 14 °) 



-a 2 /2)(p' +p' +p' +p') = - fl2 /6; (141) 



\-y ' "—|( ' f—x,+y ' f—x,—y 

for A; 3 : equivalent to Eq. (j!40P ; 
for A; 4 : equivalent to Eq. (114ip ; 



for A; 2 : 



for /c^fc y : 



for fc^fcj: 



for A; 4 .: 



- a 2 (l/2 + d/aW^+j/^+j/^) = -a 2 /Q; (142) 
ia 8 (l/2 + d/a)0/_ |4l , - pU,-») = 0; ( 143 ) 
(a 4 /2)(l/2 + d/a)(p'_ x , +y + p'- x ,- y ) = a 4 /36; (144) 



(a 4 /12)(l/2 + 2d/a - A(d/aY)(p'-x + P-x,+ y + P-x,- y ) = a*/72. (145) 

It is convenient (but not necessary) to use the normalization conditions for the probabilities. 
For the boundary sites, 

P0 + P+y + P-y + P+x,+y + P+x + P+x,-y =Po+ P+y + P-y + 1/6 = 1; (146) 

for their neighbors [sites (1,/)]) 

P- x + P'- x ,+y + P-x-y +P0+ P+x + P+y + P-y + P+x,+y + P+x-y = 

p'. x + pL Xi+y + p'. x ^y + 5/6 = 1. (147) 
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Comparing Eqs. (j!47p and (I142p . we get d = a/2, as in ID, which makes Eq. ( 1145 j) an 
identity. Then from Eqs. ffT43]) and flEHJ), 

P-«,-„ = V36 (148) 

and from Eq. (11471) . 

PU = 1/9- (149) 
Also using Eq. (fT4"g]l . from Eqs. ffT4TJ]) and (TUT]) 

P' + y=P'-y = 5/36 (150) 

and finally, from either Eq. (11391) or Eq. (11461) . 

Po = 5/9. (151) 

While 

P'- X = P- x , (152) 

P-x,+y = P-x,+y, (153) 
pi x - y = P-x-y, (154) 

the other three probabilities are different from the bulk values: 

p' +y = P+y + 1/36, (155) 
p'_ y = p- y + 1/36, (156) 
p' = po + 1/9. (157) 

Note that not just the probability of staying put changes, so unlike in ID, moves into the 
boundary are not always simply rejected. In fact, Eqs. (I155p - (I157I) can be rewritten as 

P+y = P + y+P- X ,+y, (158) 

P-y = P-y+P-x,-y, (159) 
Po = Po+P-x- (160) 

Therefore, the correct changes of the probabilities can be obtained, if the moves to the left 
(— x) from the boundary sites are rejected, but the moves in the (— x, +y) and (— x, —y) 
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directions are replaced by the moves in the +y and — y directions, respectively. In other 
words, it is as if the moves into the wall are replaced by their projections along the wall. We 
also note that, similarly to ID, the master equations for boundary sites become equivalent 
to those for bulk sites when a column of fictitious sites behind the boundary is introduced 
(open squares in Fig. |H]) and the values of the mean particle numbers at these sites are made 
equal to those of the real boundary sites across the boundary (i.e., = n (o,i))- All the 

probabilities could have been obtained from this consideration alone. It is also possible to 
obtain the 2D algorithm with a boundary as the direct product of the ID algorithm with a 
boundary and the ID algorithm in free space. 

C. Three dimensions 

We now consider the 3D case. Since we have considered two different 3D LMC algorithms 
in free space, the "direct product" algorithm with simultaneous moves in three directions 
and the unique optimal algorithm without such moves, we show how to treat a reflective 
boundary for each of these two algorithms in turn. We assume that the x = plane serves 
as the boundary. 

1. The direct product algorithm 

The extension of the direct product algorithm to the case when a boundary is present 
is straightforward and can be obtained as the direct product of the ID algorithm with a 
boundary and two ID algorithms (or one 2D algorithm) in free space. This leaves the 
probabilities of moves into bulk sites unchanged compared to the free space algorithm, but 
the "primed" probabilities for the moves into boundary sites are now 



P-s = 2/27, 



(161) 



P'± z = 5/54, 



(162) 



P—x,±y P—x,±z 



1/54, 



(163) 




,± z = 1/216, 
Po = 10/27. 



5/216, 



(164) 



(165) 



(166) 
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2. The algorithm without simultaneous moves in three directions 



In this case, the consideration is similar to 2D. Compared to the 2D case, the master 
equation has additional terms involving p' ±z , pL x ± z , and p'± y ± z - The solution of the con- 
tinuum equation is obtained from Eq. ( I137p simply by introducing the same dependence of 
the integrand on k z as its dependence on k y and integrating over k z from — oo to oo. The 
matching equation dl38[) becomes 



exp(-a 2 (^ + k 2 y + k 2 z )/6) ~ p' + p' +y e'^ a + p'_ y e^ a + p' +z e~ ik * a + p'_ 



, / -ik y a-ik z a , / -ik y a+ik z a , / ik y a-ik z a , / ik y a+ik z a 

~!~P+y,+z c "T P+y,-z c T P-y^+zt T P-y- 



e 



+(p'- x + P'-„ +y ^ a + P'- X ,- V e^ a + P'- X>+Z e~ ik ° a + ^,V fc ' B ) roS ^y ) - (167) 

First of all, there are new types of terms in the expansion with no analog in 2D, namely, 
those involving y and z (k y k z , k y k z , k y k 2 z1 k y k 2 z) k y k z , and k y k z ). Only four of the resulting 
equations are independent and they involve only four probabilities (those corresponding to 
simultaneous moves along y and z axes). Moreover, these equations are the same as they 
would be for a site in the bulk, since only moves parallel to the boundary are involved. Thus 
the four probabilities involved are the same as for bulk sites: 

P+y,+z = P'+y,-z = P'-y,+z = P-y-z = V 36 " ( 168 ) 

Equations fl 142 j) and (I145P will have the probabilities of the (—x,±z) moves added to the 
sums of the probabilities entering these equations. The normalization condition [Eq. ( I147P 
becomes 



V- x + P-x,+y + P-x-y + P-x,+z + P-x,~z + PO + P+x + P+y + P-y + P+z + P-z 

^P+x.+y + P+x,-y + P+x,+z + P+x,-z + P+y,+z + P+y,~z + P~y,+z + P-y-z = 

VL X + pL Xt+y + pL x - y + pL Xt+z + p'- x - z + 5/6 = 1.(169) 

Just as in 2D, all three equations [the modified Eqs. (I142p and (I145P and Eq. fll69p ] involve 
the same sums of "primed" probabilities, which can thus easily be eliminated, and d = a/2 
follows naturally. Equations (1 143 ft and (1 144ft remain unchanged, and there will also be 
analogous equations involving the z axis instead of the y axis, so 

p'- Xj+y = pL x - y = p'. X)+z = p'- x ^ z = 1/36. (170) 
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Then from Eq. (I169p . 

p'_ x = 1/18. (171) 

Equations (I140p and (11411) will have the probabilities of the (±?/, ±z) moves added to them 
(or subtracted in the case of the (— y, ±z) moves in the first of these equations). From these 
equations, 

p> +y =p'_ y = 1/12. (172) 

Equations obtained by matching k z and k 2 z terms are analogous to Eqs. (I140p and (I14ip . and 
from them similarly 

V' +Z = P'- Z = 1/12. (173) 
Finally, Eq. (I139p is modified to include the sum of all "primed" probabilities, and from it, 

Po = 7/18- (174) 

All other equations (modified as appropriate) and their analogs involving the z axis become 
identities once the above values are inserted in them. 



3. Relations for the probabilities 

We note that in both 3D algorithms the probabilities of all allowed diagonal moves and 
p'_ x remain unmodified compared to the bulk. All other probabilities have changed and can 
be written (again, for both 3D algorithms) as 

P+y = P+y+P-x,+y, (175) 

P-y = P-y + P-x-y, (176) 

p' +z = p +z +p- x , +z , (177) 

p'. z = p- z +p- x - z , (178) 

Po = Po+P-x- (179) 

Just as in 2D, this can be interpreted as the diagonal moves across the boundary being 
projected along the boundary, rather than rejected. Also, as in ID and 2D, the equivalence 
of the master equations for boundary and bulk sites once fictitious sites behind the boundary 
are introduced still holds, again, for both algorithms. 
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VI. ABSORBING BOUNDARIES 



Another case of practical interest is that of absorbing boundaries. In this case, particles 
reaching the boundary are "absorbed" and disappear from the system. One reason for the 
importance of this problem is that the FPT problem can be reduced to it: the cumulative 
FPT distribution is equal to the fraction of absorbed particles as a function of time. In the 
continuum diffusion problem, the absorbing boundary condition is simply 

n = (180) 

at the boundary. In LMC, since the particle number is no longer conserved, there will 
inevitably be a new type of "move" for particles in sites adjacent to the boundary, during 
which the particle simply disappears. The probability of this occurring does not enter the 
master equations explicitly (the form of these equations remains unchanged, although the 
values of the probabilities change), but the probabilities of all other moves will no longer 
sum up to unity, so the normalization conditions for boundary sites cannot be used (but 
those for bulk sites are still valid). 



A. One dimension 



The general solution of the continuum diffusion equation with an absorbing boundary at 
x = is 

POO 

n(x,t)= C{k)sm(kx)exp(-Dk 2 t)dk, (181) 
Jo 

i.e., the same as for the reflecting boundary [Eq. (1127p ], except the cosine is replaced by 
the sine. In the discrete case, the master equation for the boundary site is still given by 
Eq. (11171) . Plugging the integrand of Eq. fl 1 8 1 [) into this master equation, replacing r with 
a 2 /6D and dividing both parts by C(k) sin(kd) exp(—Dk 2 t), we get 

, ,2 2 m\ ? sin(k(d + a)) . . 

exp(-A;V/6) ~ Po ^ + Pi^o — ■ (182) 

Expanding this up to 0(k A ), we obtain 

k Qi k cl ? d ~\~ ct ( k . . k , a o, 0,0 ,q v 

1 - r — V— — = Po^o+Pi^o— 3— 1 - -^ a ( 2 d + a) + —(3a 4 + I2a i d + 8a 2 d 2 - 8ad s ) 



6 72 r r d V 6 '360 

(183) 
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It is also convenient (but not necessary) to use the condition that the sum of the probabilities 
of all moves leaving site 1 is unity (since this site is not adjacent to the boundary). Given that 
two of these moves (1 — > 1 and 1 — > 2) lead to a bulk site, their probabilities should coincide 
with the bulk probabilities ipj-^j = 2/3 and Pj+i~>j = 1/6, respectively), and therefore the 
remaining probability, pi^o, should equal 1/6. Using this in Eq. (11831) and equating the k 2 
terms, we get a quadratic equation for d that has two solutions: 

d« = a; d&> = a/2. (184) 

Equating the k° terms and using Eq. (11841) and pi^o = 1/6, we get 

p£\ = 2/3; pKo = 1/2- (185) 

As mentioned, in both cases 

pfio = = 1/6. (186) 

Both of these solutions turn the equation obtained from equating the k 4 terms into an 
identity; in fact, using that equation instead of the normalization condition would have 
given the same two solutions. Note that for the boundary site the normalization condition 
is not satisfied: the sum of the probabilities is 5/6 for the first solution and 2/3 for the 
second one. The complements of these sums to unity (1/6 and 1/3, respectively) can be 
interpreted as the probabilities of particle disappearance (or absorption by the boundary). 

Let us now consider the physical interpretation of the two solutions. Note first of all that 
the first solution corresponds to the absorption probability of 1/6, which is the same as the 
probability of the move to the left from a bulk site, and the probabilities of all allowed moves 
are unchanged compared to the bulk. Therefore, when using the first algorithm, we can treat 
the boundary site as a bulk site, except whenever a move to the left is attempted, the particle 
is deleted instead. There is no such simple interpretation for the second solution. Just as for 
the reflecting boundary, we can also analyze the algorithms by looking at the corresponding 
master equation. If the first set of parameters is used in the master equation (I117p . the 
resulting equation coincides with that for bulk sites [Eq. (TT7]) ]. if a fictitious site (numbered 
— 1) is introduced and the mean particle number at that site is made identically equal to 
zero. The physical interpretation of this is clear: if d = a, then site —1 is placed right at the 
absorbing boundary (Fig. [9]), where the particle concentration in the continuum problem is 
indeed exactly zero. On the other hand, if the second set of parameters is used in the master 
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FIG. 9: One of the two variants of treatment of an absorbing boundary in the ID optimal algorithm. 
The inaccessible region is shaded. The site closest to the boundary is the boundary site (filled 
square); all other sites are bulk sites (filled circles). In this variant, the placement of the boundary 
site is at the distance d = a from the wall. The open square is a fictitious site at the wall useful 
in the analysis of the algorithm, as described in the text. The move left from the boundary site 
(dotted line) leads to the disappearance of the particle. Only the probabilities of the moves to the 
boundary site (thick arrows) are allowed to change in the analysis, but turn out to be the same as 
in the bulk for the optimal algorithm. Choosing d = a/2, as in Fig. is also possible, but leads to 
different probabilities of the moves. 



equation, the resulting equation will coincide with that for bulk sites, if n_i = — uq. Since 
for d — a/2 the wall is halfway between sites —1 and 0, we get n = at the wall simply by 
interpolating between no and Of course, linear interpolation can be done for any d; for 
the interpolated value at x = to be zero, 

d-d / 10 r,N 

n_i = — no- (187) 

d 

However, it is only at a = d/2 and a = d that at least the cubic correction to the interpolated 
value vanishes and thus the interpolation is more precise (in the latter case, of course, no 
interpolation is involved, as site —1 is located right at the boundary). 

Note that if we are content with the linear accuracy of interpolation, we can use Eq. I I 18 7ft 
in Eq. (JT7J and by comparison with Eq. (j!17p obtain 

a — d 

Po^o = Pj^j ~ Pj-i-ij • (188) 

ct 

Note that in the ordinary algorithm, Pj^j = and so po->o is non-negative and thus physically 
meaningful only for a = d. Thus there is only one way of implementing the absorbing 
boundary condition within the framework of the ordinary algorithm. Introducing a nonzero 
probability of staying put broadens the set of allowed values of d and thus gives more 
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freedom. In particular, for Pj-+j = 2/3 and pj-i^j = 1/6, as in the optimal algorithm, any 
d > a/5 is possible. Even when optimality is required, we still have two choices, d = a/2 
and d = a, more than for the ordinary algorithm. This is yet another advantage of using 
the optimal algorithm, or nonzero Pj->j in general. 

Given the choice between the two optimal algorithms (with d = a and d = a/2), which 
one should be preferred? In principle, the first choice has certain advantages. First, it 
follows from the previous discussion that the boundary condition n = is explicit in the 
first algorithm, but is obtained indirectly via interpolation in the second. It is not obvi- 
ous, however, that this improves the accuracy of the algorithm overall, as in both cases 
the continuum solution satisfies the respective master equations to the same order [0(/c 4 )]. 
Second, recalling the correspondence between diffusion with absorbing boundaries and the 
FPT problem, we note that it is the variant with d = a that translates directly into the 
discrete FPT problem as considered in Sec. Ill CI in this case the rate of reaching a site is 
equal to the absorption rate when that site is replaced by an absorbing wall. This means 
that this variant of the algorithm will produce two correct moments of the absorption time 
distribution; it can be checked that this is not the case for the other variant. However, there 
can be other situations when the d = a/2 variant is preferable. 

B. Two dimensions 

As in the case of a reflecting boundary, we assume that x = is the boundary and x > 
is the accessible region. We can derive the parameters of the optimal algorithm using the 
same approach as in the reflecting boundary case, by writing down the continuum solution, 
which in this case is 

/*oo /*oo 

n(x,y,t) — / / C(k x , k y ) sm(k x x) exp[ik y y — D(kl + ky)t]dk x dk y , (189) 

J k x =Q J k y =—oo 

plugging the integrand into the master equation (JT36J) and demanding that the resulting 
equality is satisfied as accurately as possible. A much shorter route is obtaining the optimal 
algorithm as the direct product of the ID algorithm with an absorbing boundary and another 
one in free space. Since there are two variants of the optimal ID algorithm (with d = a and 
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with d = a/2), we end up with two variants of the 2D algorithm: variant 1 - 

d (1) = a, (190) 

Po (1) = 4/9, (191) 

V { -l=V± y = 1/9, (192) 

p'W ±y = 1/36, (193) 

and variant 2 — 

= a/2, (194) 

Po (1) = 1/3, (195) 

= 1/9, (196) 

p'S = 1/12, (197) 

p'W ±y = 1/36. (198) 

The particle disappearance probabilities are 1/6 and 1/3, respectively, as in ID. Note that 
in the first variant all probabilities are the same as the bulk probabilities. Therefore, as 
in ID, the boundary sites can be treated as the bulk sites, except moves into the walls are 
replaced by particle disappearance. This is not the case for the second variant. Also, it is 
still true, as in ID, that the master equations for the boundary sites coincide with those for 
the bulk sites with n_i j = for the first variant and n_i j = — n j for the second variant. 

The disadvantage of the direct product approach is that, strictly speaking, we do not prove 
that the resulting sets of parameters are the only possible ones. A much longer consideration 
using the first method described above shows that the two variants of the algorithm we have 
obtained are indeed the only sets of parameters turning the master equation into an identity 
up to 0{k A ) when the continuum solution [Eq. (11891) ] is inserted into it. 

C. Three dimensions 

In 3D, both in free space and with a reflecting boundary, we considered two algorithms: 
the "direct product" algorithm and the algorithm with moves in one and two directions only. 
We now extend these considerations to the case of an absorbing boundary at x = 0. 
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1. The direct product algorithm 



The direct product algorithm is obtained as the direct product of the 2D algorithm with 
an absorbing boundary and a ID algorithm in free space. Since there are two variants of 
the former, we end up with two variants of the 3D algorithm. In the first variant, 

d {1) = a, (199) 

Po (1) = 8/27, (200) 

P*2 = PS = PS = 2/27, (201) 

V% y = pl X l ±z = p'l^ ±z = 1/54, (202) 

P-S±v,±* = 1/216; (203) 



in the second variant, 



d^ = a/2, (204) 

P? } = 2/9, (205) 

P-2 = 2/27, (206) 

p'i 2 J = pTJ = 1/18, (207) 

p'® ±y =P% z = 1/54, (208) 



pSL = 1/72, (209) 



P-i± y ,±z = 1/216. (210) 
2. The algorithm without simultaneous moves in three directions 

To quickly derive this algorithm, we use an approach based on the correspondence between 
the master equations for the boundary and bulk sites, as described below. 

First, we have mentioned when discussing the ID algorithms in Sec. IVI Al that for the 
variant with d = a these master equations coincide when n_i is put equal to zero in the 
bulk equation. Extending this to 3D, we replace ri_i^, m with zero for all / and m in the 3D 
bulk master equation for the evolution of n j 0tmo and compare to the master equation for a 
boundary site. The result is that all probabilities remain as in the bulk [i.e., are given by 
Eqs. (1109]) — (lllip ], except that moves involving the +x direction are, of course, impossible; 
also, the complement of the sum of the probabilities to unity (which equals 1/6, as in ID 
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and 2D) is the probability of the particle disappearance. The same interpretation as in ID 
and 2D is still valid: the boundary sites can be treated as bulk sites, but all attempts to 
move into the boundary are replaced by particle disappearance. 

On the other hand, for the variant with d = a/2 the ID master equations for the boundary 
and bulk sites coincide when n_i = —n Q in the latter. Again, extending this to 3D, we obtain 
the following probabilities: 

Po=Po-P +x = 5/18, (211) 
p'„ x =p- x = 1/18, (212) 

P±y = P±z = P±y - P+x,±y = P±z ~ P+x,±z = 1/36, (213) 
P'-x,±y = P-x,±z = P±y,±z = P~x,±y = P-x,±z = P±y,±z = 1/36. (214) 

The sum of all these "primed" probabilities is 2/3, and so the disappearance probability is 
1/3, as in ID and 2D. 



VII. DISCUSSION 



In this paper, we have obtained the sets of parameters (transition probabilities and the 
time step) of LMC algorithms for particle diffusion problems that optimize the accuracy of 
these algorithms. The problem was solved by demanding that the solutions of the master 
equation of the algorithm approximate those of the continuum diffusion equation as accu- 
rately as possible. The matching between the continuum equations and the discrete master 
equations was done using the general solution of the continuum equation written as a Fourier 
expansion in space where each component decays exponentially in time. The general solu- 
tion of the master equation written in a similar form can also be used, in which case the 
goal is to make the decay rates of the Fourier modes as close as possible in the continuum 
and discrete equations. However, having the solution of the master equation is not neces- 
sary: for example, we did not use it in Sees. |V] and IVIt where we considered the treatment 
of boundaries. In that case, the solution of the continuum equation was simply inserted 
into the master equation and the parameters of the latter adjusted to satisfy the resulting 
equality as accurately as possible. Having the solution of the master equation, while not 
necessary, allowed the illustrative comparison between the dispersion relations presented in 
Fig.d 
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Since LMC simulations are generally carried out in order to study long length and long 
time scale diffusion processes, all expressions are expanded in the wavenumber k and terms 
of the lowest order are matched. Matching the dispersion relations up to 0(k 2 ) guarantees 
that the average velocity (equal to zero in the unbiased case) and the diffusion constant are 
reproduced correctly by the algorithm. In other words, both the mean displacement (the first 
moment of the distribution for particles starting at a particular site; again, equal to zero) and 
the mean-square displacement (the second moment) remain correct at all times. Moreover, 
for all higher-order moments, the leading term at large times has the correct coefficient. 
However, these higher moments may be incorrect for shorter times. The algorithms can 
suffer from some artifacts, as discussed in Sec. [Ill for the ID case. Going to 0(k 4 ), which 
is possible with only first neighbor moves in ID and when both first neighbor and second 
neighbor (diagonal) moves are included in 2D and 3D, removes the artifacts and makes the 
fourth moment, as well as the subleading terms of even higher moments, correct, too, which 
we have shown for the ID case, but expect to generalize to 2D and 3D as well. We refer to 
the resulting algorithms as optimal, in the sense that they achieve the best accuracy given 
the set of moves, although they are not optimal in terms of the simulation speed if accuracy 
is not a concern. The latter fact is due, in particular, to a distinctive feature of the optimal 
algorithms: a nonzero probability for a particle not to move during a particular time step, 
or a waiting time. 

We have also shown for the ID case that for a particle starting at a particular point 
between two boundaries, the optimal algorithm reproduces correctly the first two moments 
of the distribution of the times of first passage to the boundaries, if both the initial point and 
the boundaries coincide with lattice sites. In fact, the optimal algorithm can be obtained 
based on the requirement that the first two moments of the first-passage time are reproduced 
correctly or on the requirement that the first four moments of the particle distribution 
are correct, instead of matching the dispersion relations. We have also studied the full 
distributions of the first-passage times comparing them to the exact continuum result. We 
have shown that the optimal algorithm converges much faster to the exact result than the 
ordinary algorithm without a waiting time as the mesh step decreases, so if a particular 
accuracy is required, a much coarser mesh can be used with the optimal algorithm, which 
provides a speedup that in the end more than compensates for the loss of computational 
speed due to the waiting time. 
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One fact worth mentioning is an analogy between the optimal LMC algorithms and 



Lattice Boltzmann (LB) algorithms for mesoscopic fluid simulations 24|, |25|. LB equations 
can be thought of as master equations for particles residing at the sites of a lattice and 
moving at each time step to a predetermined set of nearby sites, as in LMC. However, a 
crucial difference is that the probabilities of particular moves are not fixed. At each site, 
the mean numbers of particles with each particular velocity are considered separately and 
evolve individually according to the LB equations. The velocity of a particle determines 
how that particle moves at the next time step, and since the set of moves is discrete, the 
set of possible velocities in the algorithm is discrete as well. In different varieties of the 
LB approach different sets of discrete velocities (and thus of possible moves) are used. The 
standard notation used to distinguish these varieties is DmQn, where m is the dimensionality 
of the space and n is the number of possible velocity vectors. One question relevant to 
the design of LB algorithms is finding the best approximation of the equilibrium Maxwell 
velocity distribution using a particular set of allowed discrete velocity vectors. It turns out 
that in the best approximation, the mean numbers of particles with particular velocities are 
proportional to the probabilities of the corresponding moves in the optimal LMC algorithm 



with the same set of moves. For the D2Q9 LB algorithm 



261 ]. this will be our 2D optimal 



LMC; for the D3Q19 algorithm 25j, our 3D LMC with moves along one and two directions; 



for the D3Q27 algorithm 



26l | . our 3D "direct product" algorithm. This correspondence is 



not coincidental. The Maxwell distribution is a Gaussian distribution, as is the distribution 
of positions of particles diffusing in continuum space at a given time after they started at 
the same point. Given this, it is easy to see that the problem of finding the distribution of 
discrete velocities approximating the Maxwell distribution is mathematically equivalent to 
finding the set of probabilities of moves such that after a single step starting from the same 
site the resulting particle distribution is as close to Gaussian as possible. Of course, even 
though there is a mathematical equivalence, the physical meaning of the problems is very 
different. 

Another interesting result of this paper concerns the treatment of impenetrable reflecting 
boundaries. Often in MC algorithms, if an attempted move is forbidden, e.g., because a 
particle would overlap with an obstacle, it is simply rejected. This is a correct choice in 
algorithms that only intend to reproduce equilibrium properties, as it preserves detailed 
balance, which is all that matters in that case. However, as we have shown, this is not the 
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best choice in a dynamical algorithm in a two- or higher-dimensional space. In that case, 
the best accuracy is achieved by replacing forbidden moves across the boundary with their 
projections along that boundary. This result makes physical sense: if, for example, in 2D all 
boundaries are orthogonal to the x direction and are infinitely long in the y direction, the 
diffusion in the y direction should not be affected, which is only possible if the simultaneous 
moves in the x and y directions are projected along the y direction. 

We have also considered the case of absorbing boundaries. It is interesting that for each 
of the optimal free space algorithms, there are two ways of treating the boundaries parallel 
to one of the lattice axes. The more obvious one is to destroy the particle whenever it moves 
into the wall, but leave all other moves unmodified. This is the best approximation for a 
wall that is at the distance equal to the mesh step from the last row of sites next to the 
boundary (d = a using the notation of the paper). However, there is also a second variant, 
with rules that do not have a simple interpretation, which corresponds to the last row of 
sites being at the distance from the boundary of half the mesh step (d = a/2). This second 
variant (as well as intermediate, non-optimal choices) is only possible in algorithms with a 
waiting time (even if optimality is not required). 

Note that we have only treated the case of infinite, flat boundaries. In real situations 
of interest, one deals with finite and/or curved boundaries (finite obstacles, tortuous pores, 
etc.) In the simplest case, the walls consist of flat pieces, each of which is orthogonal to 
one of the axes, that are joined at corners, like those in Fig. [10J An example would be a 
rectangular obstacle in 2D with the sides along the axes. Sites adjacent to the walls that 
are far from the corners can be treated as boundary sites next to infinite walls (this assumes 
that it is possible to choose the lattice so that the distance d from the boundary sites to 
the boundary is as required by the algorithm everywhere). However, we still need a way to 
treat sites next to the corners. We note, though, that we have assumed throughout that the 
particle concentration varies smoothly on the length scale of the mesh step, and this implies 
that any features of the obstacles and the walls are likewise smooth on this length scale, so 
corners should be rare. It does not matter much then how exactly these rare situations are 
treated; we will consider these complications in a separate paper in the future. However, we 
can offer some intuitive guidance in some of these situations. For example, it is intuitively 
clear how to deal with absorbing boundaries in the algorithms with d = a [Fig. [TTJJ, (a) and 
(b)]: in these algorithms, all moves into the wall lead to absorption and all other moves are 
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unchanged, which generalizes straightforwardly to corners. It is also clear how to handle the 
situation with d = a/2 (either absorbing or reflecting walls), when, e.g., in 2D the particle 
is allowed to move in one quadrant and the other three are blocked [Fig. [TU] (c)]. In this 
case, fictitious sites are introduced inside the forbidden area symmetrically with respect to 
the boundary to the row of sites nearest to it [Fig. [TUJ (c)], and the mean particle numbers 
at these sites are assumed equal to those at sites nearest to them across the boundary (with 
the minus sign in the case of absorption). The generalization to the case of a corner is 
straightforward, by assuming that the fictitious site closest to the corner takes on the value 
of the real site closest to the corner (possibly with the minus sign). However, in the otherwise 
analogous case when three quadrants are free and one is blocked [Fig. []Jj] (d)], there is some 
ambiguity, as it is not clear which of the three real sites closest to the corner should give its 
value to the fictitious site nearest the corner. This translates into an uncertainty about the 
probabilities of the diagonal moves "passing through the corner" shown in the figure. 

Curved boundaries or those not orthogonal to the axes are even more complicated. Ap- 
proximating such boundaries by piecewise flat boundaries with each piece orthogonal to one 
of the axes, in such a way that those sites inside the original boundary are still inside and 
those outside are still outside, would in effect produce "jagged" boundaries with many cor- 
ners adjacent to each other that need special treatment. A more accurate approach would 
be approximating such boundaries with flat pieces oriented arbitrarily. However, this also 
requires a separate consideration. The rules we have derived cannot be extended straight- 
forwardly to such "skewed" boundaries, even when they are perfectly flat: for instance, 
projecting forbidden moves along the boundary is impossible, as such a projection would 
not coincide with any of the allowed moves. Solving these problems will be a subject of 
future work. Considering finite hard particles of an arbitrary shape is a more complicated 
problem, since particle rotation needs to be taken into account (spherical particles are an 
exception). An even more complicated problem is that of several diffusing hard-core par- 
ticles that can collide with each other. Such particles would move independently between 
collisions (according to the rules derived in this paper, if rotation can be ignored), but rules 
for treating collisions still need to be worked out. 

The preceding discussion may sound as if there is still some way to go before the al- 
gorithms described here become practically useful. This is not entirely true, however. As 
mentioned in Sec. Ill Al our ID optimal algorithm is the zero-field limit of a previously de- 
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FIG. 10: 2D examples of different configurations of walls forming corners and different variants of 
meshes for optimal LMC algorithms. Forbidden areas are shaded; real sites are black circles and 
fictitious sites introduced to facilitate the analysis of the algorithms are white circles. Cases (a), 
(b) and (c) can be analyzed, as described in the text, but in the case (d) there are complications 
that physically correspond to the uncertainty about the probabilities of the diagonal transitions 
"across the corner" indicated by the arrow. 



rived algorithm for biased diffusion 19j. We have recently applied the latter algorithm to a 



ID model of polymer translocation through a nanopore [27| and showed that this algorithm 
gives very different widths of the translocation time distributions compared to the ordinary 
algorithm without waiting times, thus demonstrating the importance of introducing waiting 



times [28[. It is true, though, that the utility of the new algorithms will be further enhanced 
when the issues described above are taken care of, as well as when the approach is extended 
to biased diffusion in 2D and 3D as well. In general, we expect the new algorithms to be 
useful in those problems where the transient behavior of the system at short times is im- 
portant. This includes various first-passage problems (including the polymer translocation 
problem mentioned above), as well as adsorption, aggregation jl], anomalous diffusion 
and reaction-diffusion {3]] problems. One example of a problem dealing with transient effects 
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where ou r appro ach may be useful is studying drug release profiles from disordered porous 



matrices 



151 ] . On the other hand, if only the long-time behavior is of interest and short- 



time transient dynamics does not affect it, then the ordinary algorithms with jumps at every 
time step may be just as good and may offer computational time savings compared to the 
optimal algorithms. 

In those cases where discrete-time LMC can be used, there is also an option to use MC 
algorithms with continuous space (off-lattice MC), continuous time (often referred to as 
kinetic MC £ , S 31]) or both. Of course, we do not claim that LMC is always the best 
option; rather, the goal of this paper is to suggest that in those cases where the choice 
in favor of LMC is made, one might consider improving the accuracy by using the LMC 
modifications described here. We do note, however, that since the choice between LMC and 
continuous MC is often driven by the tradeoff between higher speed of the former and higher 
accuracy of the latter, improving the accuracy of LMC may actually make LMC preferable 
in some situations where this would not be the case otherwise. 
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Appendix: a derivation of the long-time dependence of the first-passage rate for 
the LMC algorithms 

The ID first-passage problem for LMC algorithms, as formulated in Sec. Ill C\ can be 
solved by considering the master equations ( ITT]) for sites from —N + 1 to iV — 1 and fixing 

n ±N = 0. (A.l) 

The solution of this system will give the mean numbers of particles at every site that have not 
yet reached the walls at sites ±N. This solution can be obtained from the general solution 
for the infinite lattice [Eq. (ITS]) ] by picking those modes that satisfy the boundary conditions 
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[Eq. flA~Tj) ]. which gives k m = vr(2m + l)/2Na with m = 0, . . . , N — 1 and C(-A;) = C(fc); 



thus the solution is 

N-l 

nj (t) = ^C m cos(7r(2m + l)j/2N) exp[-a d (k = vr(2m + l)/2Na)t]. (A.2) 

m=0 

The mean number of particles reaching the walls at time step occurring at time t, v(t), is 
proportional to the number of particles at adjacent sites at the previous time step, i.e., 

v{t) = Pj-t^jUN^t - r) + p j+1 ^jn- N+ i{t - r). (A. 3) 

This number per unit time (i.e., the rate) is 

r{t) = = - T %_ 1 t - r , A.4 
r cr 

where we have used Eq. (fTTj) and the fact that, according to Eq. (IA.2I) or simply from 
symmetry considerations, n_N+i = n N-i- 

For the optimal algorithm, ad(fc) is a monotonically increasing function, thus in the large 
t limit, only the first term in the sum survives, which gives 

rij(t) ~ C cos{irj/2N) exp[-a d (k = ic/2Na)t]. (A.5) 



Then from Eq. (1A.4I) 



r(t) ~ — sin(7r/2iV) exp[a d (fc = h/2No)t] exp[-a d (k = n/2Na)t}. (A.6) 
a 2 

Therefore the decay rate is 

QDN 2 

op t = a d {k = n/2Na) = — ln[2/3 + (1/3) cos(tt/2A0] , (A.7) 

tr 

where we have used a = b/N and Eq. (f!9l) with p and r given by Eqs. (I2"6"l) - (128j) . The 
prefactor is 



2DC ( 



7o P t = sm(n/2N) exp[a d (k = n/2Na)r]. (A.8) 



a 

To find C , we need to expand the initial condition, rij(0) = 5j )0 in the Fourier series, 

N-l 

S jfi = ^ C m cos(7r(2m + l)]/2N). (A.9) 

m=0 

Since 

N 

cos^^ + l)]/2N) cos(vr(2m 2 + l)]/2N) = NS muma , (A.IO) 

j=—N 
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A' 



we get 

C ™ = Jj E ^ocos(vr(2m + l)j/2N) = i 

j=-JV 

for all m, including 0, and therefore 
2D 



7o P t 



sin(vr/2iV) exp[a d (k = ti/2No)t] 



2DNsm(7r/2N) 



(All) 



(A.12) 



A^a 2 ' y J 6 2 [2/3 + (l/3)cos(7r/2AT)]' 

For the ordinary algorithm, the situation is slightly more complicated, because a d (k) is 
no longer monotonic; as a result, |ad(^iV-i)| = a d(^o) an d two terms of equal magnitude 



survive at large time: 



rijit) ~ C cos ( — ) exp 



irj 



+Cn-i cos 
The first-passage rate is 
r(t) 



a d ( k 
n(2N - l)j 



7T \ 



2Na, 



2N 



exp 



-a d k 



2D ( 7T \ 

iVV Sin l2ivJ 



exp 



-a d (A; 



7T 



+ (-l) 7V+1 exp 



-ot d k 



2Na 
tx{2N - 1) 

2iV^ 



tt(2N- 1] 
2iV^ 

(t-r) 
(t-r 



(A.13) 



(A.14) 



where we have used Eq. (lA.lljl for C m (that is still valid for the ordinary algorithm). Ac- 
cording to Eq. ( |T9l) with p and r given by Eqs. (JIJ) (J3j) , 



ra d (k 



ra d k 



2Na 
tt(2A - 1 



) 



2Na 



— In cos ( 
-In 



cos 



\2NJ 



ra d [k 



71 \ 



2Na, 



ITX, 



so Eq. ( 1A.14I) becomes 
2D 



r(t) 



. ( IT \ 

sin 

No? V2iV/ 



exp 



-a d Ik 



IT 

2iV^ 



(t 



x 



[1 + 



(A.15) 
(A.16) 

(A.17) 



where M is the number of the time step (t = Mr). The expression in the rightmost square 
brackets alternates between and 2, which represents a familiar feature of the ordinary 
algorithm: no particles reach the boundaries at even steps for odd A" and at odd steps for 
even N. Since we have decided to ignore these oscillations, we simply replace this expression 
with its average value of unity. Then we immediately obtain for the decay rate, 

2DN 2 



ord 



a d (k 



7T \ 



) 



2Na 

Tord 



111 cos iw) 



b 2 

2DN ( 7T \ 



(A.18) 
(A.19) 
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